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1.  Introduction  and  Summary 

The  matching  or  overlaying  of  images  from  similar  or  dissimilar  sensors 
is  termed  "registration".  Matching  of  images  with  other  images  including  sym¬ 
bolic  images  such  as  maps  is  a  general  and  important  problem  in  a  wide  variety 
of  image  analysis  tasks,  in  particular  in  those  aimed  at  extracting  mapping, 
charting  and  geodetic  data.  The  registration  of  digital  and  digitized  images 
and  digital  maps  is  necessary  in  planned  systems  for  creating  and  updating 
digital  cartographic  data  bases  such  as  the  Digital  Landmass  System  (DLMS) 
database. 

Registration  is  distinguished  from  rectification  which  refers  to  the  at¬ 
tempted  correction  of  geometric  distortions  introduced  during  the  acquisition 
of  the  imagery  by  the  sensing  mechanism  and  a  variety  of  other  sources  which 
are  categorized  and  briefly  described  in  section  2. 

As  noted  in  section  2,  not  all  distortions  are  correctable  without  regis¬ 
tering  the  imagery  with  a  database.  This  interdependence  between  rectification 
and  registration  together  with  the  different  resolutions,  sensing  geometries 
and  different  appearances  of  contents  of  imagery  acquired  by  different  sensors 
makes  registration  a  difficult  problem.  In  general,  however,  most  imagery  of 
interest  will  depict  distinguishing  structural  characteristics  such  as  lines, 
shapes,  or  tones,  the  aggregate  of  which  allows  the  identification  and  analysis 
of  the  scene.  While  the  individual  structural  features  may  be  interpreted 
ambiguously  the  global  geometric  arrangement  often  admits  an  unambiguous  inter¬ 
pretation  for  the  scene  and  its  various  regions. 

Individual  points  in  an  image  field  cannot  always  be  uniquely  interpreted. 
Points  which  are  not  uniquely  identifiable  by  local  structure  and  global  geo¬ 
metry  are  called  ambiguous  points.  A  pass  point  is  a  point  in  an  image  iden¬ 
tifiable  by  the  structure  of  its  local  neighborhood  and  its  geometric  relation- 
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ship  to  other  pass  points  and  identifiable  structures.  Examples  are  forks  in 
streams,  mountain  peaks,  or  road  intersections.  Points  in  wheat  fields  or 
lakes  are  not  usually  pass  points.  Pass  points  whose  absolute  location,  i.e., 
latitude,  longitude  and  elevation,  on  the  globe  is  known  are  called  ground 
control  points. 

Fortunately  there  are  usually  enough  pass  points  evident  in  reasonably 
rich  imagery  to  permit  global  matching  of  images  even  without  ground  control  or 
platform  attitude  information.  Automatic  procedures  for  global  matching  based 
on  pass  point  determination  must  however  confront  the  fact  that  pass  points  are 
very  sparsely  distributed  and  the  bulk  of  the  image  points  are  ambiguous  points. 

This  document  reports  on  our  study  of  the  problem  of  automatic  registration 
of  images  with  maps  and  images  from  dissimilar  sensors.  The  registration  prob¬ 
lem  is  divided  into  four  subproblems:  (a)  detecting  and  extracting  appropriate 
features  for  pass  point  determination;  (b)  determination  of  an  approximate 
global  registration;  (c)  determination  of  disparity  for  a  subset  of  points  in 
the  images;  (d)  determination  of  the  global  nonlinear  transformation  for  the 
entire  image. 

For  all  methods  other  than  the  standard  correlation  of  gray-scale  imagery, 
the  inexpensive  extraction  of  good  features  is  a  primary  concern.  The  cost  of 
feature  extraction  must  be  balanced  against  the  potential  use  of  the  features. 
Region  features,  while  often  costly  to  obtain,  can  be  vital  in  later  phases  of 
image  interpretation.  On  the  other  hand,  isolated  interest  points,  such  as 
points  of  high  curvature,  can  be  inexpensive  to  obtain  but  may  serve  no  purpose 
other  than  to  determine  registration. 

Following  a  brief  discussion  on  rectification  in  section  2,  section  3  con¬ 
siders  the  detection  and  extraction  of  edge  features,  point  features  and  region 
features.  Edge  detectors  based  on  gradient  and  Laplacian  operators  and  Hough 


transforms  are  described.  A  simple  linking  method  which  produces  reasonably 
long  edge  segments  and  produces  good  results  is  also  presented.  Point  feature 
detectors  include  two  point  operators  and  detectors  for  intersections  and  high 
curvature  points.  Region  features  based  on  invariant  moments  and  region  seg¬ 
mentation  are  described.  This  section  also  presents  results  of  applying  some 
of  the  feature  detectors  to  a  panchromatic  image  and  infrared  image  of  the  same 
scene. 

Section  4  is  devoted  to  approaches  to  registration  and  includes  a  brief 
survey  of  some  existing  techniques  for  matching.  Section  4.1  describes  the 
basic  correlation  method,  a  sequential  correlation  method,  a  hierarchical  cor¬ 
relation  method,  presents  example  of  correlating  feature  vectors  where  the  fea¬ 
tures  are  edges,  and  another  example  using  abstract  triangles  formed  by  con¬ 
necting  three  successive  high  curvature  points.  Section  4.2  presents  the  basic 
L.N.K.  registration  procedure  and  examples  of  its  application  to  images.  In¬ 
cluded  is  a  brief  report  on  an  experiment  with  an  example  involving  scale  change. 
The  experiment  demonstrates  the  utility  of  the  proposed  technique  for  regions 
where  cultural  activity  creates  features  such  as  straight  edges  or  networks 
of  lineals.  The  results  reported  are  typical  of  many  similar  experiments. 

Section  4.2.3  presents  an  extension  of  the  LNK  registration  method  to  3- 
dimensional  images.  Good  approximate  rotation,  scale  and  translation  (RS&T) 
transformations  were  obtained  automatically  for  photo/map  pairs  even  when  there 
was  some  relief  in  the  terrain.  However  there  are  many  cases,  such  as  in  low 
altitude  aerial  imaging,  where  RS&T  transformations  are  inadequate.  In  such 
cases  projective  transformations  must  be  used  and  this  section  shows  how  the 
LNK  registration  technique  can  be  extended  to  3-d  under  certain  constraints. 

All  versions  of  the  L.N.K.  registration  technique  involve  clustering  in 
parameter  space.  Section  4.2.4  discusses  two  techniques  -  hierarchical  cluster- 
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ing  and  variable  resolution  clustering  -  which  L.N.K  has  used  in  its  registra¬ 
tion  procedure. 

Symbolic  matching  of  images  consists  of  reducing  two  images  to  a  pair  of 
abstract  structures  and  comparing  these  structures.  Section  4.3  considers 
methods  which  take  the  regions  in  a  segmented  image  as  components  of  the  struc¬ 
ture  describing  the  image.  Most  of  the  available  region  matching  algorithms 
do  not  provide  the  registration  transformation  but  only  a  correspondence  be¬ 
tween  some  subset  of  the  set  of  regions  in  the  first  image  and  a  subset  of  the 
set  of  regions  in  the  second  image.  Section  4.3.1  describes  region  image  match¬ 
ing  using  similarity  of  region  features  and  section  4.3.2  considers  an  approach 
in  which  graph  theoretic  techniques  are  used  to  match  region  segmented  images 
represented  as  graphs.  Each  node  of  the  graph  corresponds  to  a  region  and  an 
edge  in  the  graph  represents  adjacency  between  two  regions.  In  this  Region 
Adjacency  Graph  (RAG)  each  node  also  contains  a  set  of  descriptors  such  as  shape, 
size,  and  average  brightness,  about  the  region  represented. 

Section  4.4  presents  a  comparison  of  registration  procedures.  Because  of 
the  unavailability  of  relevant  imagery,  the  scope  of  the  contract  and  the  less 
than  precise  descriptions  available  in  the  publications  describing  various  meth¬ 
ods,  it  was  feasible  to  present  only  a  brief  qualitative  comparison. 

Once  a  global  registration  has  been  determined  it  is  usually  necessary  to 
modify  the  transformation  to  account  for  local  distortions.  The  disparity  be¬ 
tween  two  registered  images  is  the  small  local  differences  in  corresponding 
pixel  locations  caused  by  these  distortions.  Section  4.5  discusses  the  deter¬ 
mination  of  image  disparity.  Section  5  presents  our  conclusions  and  recommenda¬ 
tions  for  further  investigation. 


2.  Rectification 


All  remote  sensing  imagery  have  geometric  distortions.  These  distortions 
can  be  grouped  into  six  categories: 

(1)  distortions  caused  by  the  topography, 

(2)  distortions  caused  by  the  sensing  mechanism, 

(3)  distortions  caused  by  the  recording  mechanism, 

(4)  distortions  caused  by  a  non-ideal  flight  path, 

(5)  distortions  caused  by  perspective, 

and  (6)  distortion  caused  by  movement  of  the  earth. 

All  these  distortions  would  make  registration  more  difficult.  Thus,  all  cor¬ 
rectable  distortions  should  be  removed  before  applying  any  registration  tech¬ 
niques. 

Distortions  due  to  topography  include  changes  in  position  caused  by  the 
terrain's  altitude  and  distortions  caused  by  the  earth's  curvature.  In  princi¬ 
ple,  the  distortions  due  to  altitude  can  be  eliminated  if  there  is  an  altitude 
database  available  for  that  region.  The  only  problem  with  using  it,  is  that  the 
image  must  be  registered  first.  Because  of  this  conflict,  registration  tech¬ 
niques  should  be  flexible  enough  to  allow  for  topographical  distortions. 

If  the  distortions  caused  by  the  sensing  and  recording  mechanism  can  be 
modelled,  then  it  may  be  possible  to  completely  remove  the  distortion.  In 
most  cases,  some  of  the  parameters  of  the  model  are  known 
approximately.  In  this  case  only  some  of  the  distortion  can  be  removed. 

For  example,  tangential  distortion  of  line  scanner  imagery,  in  which  the 
scale  varies  across  the  flight  path,  causes  lineal  features  across  the  flight 
path  to  have  an  "S"  shape.  In  order  to  precisely  determine  the  necessary  cor¬ 
rection  for  a  particular  pixel,  the  height  of  the  sensor  above  that  pixel  needs 
to  be  known.  As  mentioned  before,  even  with  an  available  altitude  database. 


the  image  must  be  registered  before  it  can  be  used.  For  the  above  example, 
the  tangential  distortion  is  removed  by  using  an  approximate  height  above  the 
terrain  which  is  usually  the  average  altitude  of  the  aircraft  carrying  the 
sensor. 

Distortions  due  to  non-ideal  flight  path  include  those  due  to  any  deviation 
from  the  optimal  orientation  of  the  aircraft  and  distortions  due  to  movement  of 
the  aircraft  during  image  taking. 

Distortions  due  to  non-optimal  orientation-  share  the  problem  that  their  cor¬ 
rections  need  the  terrain's  altitude.  Correction  requires  the  parameters  of  the 
flight  path  to  be  known. 

Distortions  due  to  perspective  occur,  in  addition,  if  the  orientation  of 
the  sensor  is  not  vertical  or  the  f ield-of-view  of  the  sensor  is  large.  Precise 
correction,  again,  requires  the  terrain's  altitude. 

Finally,  movement  of  the  earth  can  be  a  problem  in  satellite  scanner  imagery 
since  it  causes  the  image  to  have  a  shape  like  a  parallelogram. 


3.  Feature  Selection 


3. 1  Feature  Detection 

The  role  of  feature  selection  in  registration  is  an  important  one.  Each 
registration  method  requires  a  certain  image  representation,  which  contains  a 
set  of  features  representing  the  image.  This  in  turn,  means  a  set  of  feature 
detectors  to  detect  those  features  is  needed.  No  registration  technique  can 
be  better  than  its  corresponding  feature  detectors.  If  the  method  requires 
edges  and  the  edge  detector  worked  poorly,  then  the  registration  will  not  be 
good.  Hence,  what  features  can  be  detected  cheaply  and  reliably  will  greatly 
influence  the  choice  of  a  registration  method. 

This  section  is  divided  into  three  parts,  edge  features,  point  features, 
and  region  features.  In  order  to  show  the  utility  of  some  of  the  features 
presented,  the  feature  detectors  were  applied  to  a  panchromatic  image  and  in¬ 
frared  image  of  the  same  scene.  The  results  were  then  correlated  and  examined. 
The  correlation  techniques  used  were  basic  ones  and  are  later  discussed  in 
section  4.1.1. 
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3.2  Edge  Detection 


There  is  evidence  that  registration  of  two  images  should  be  based  on 
the  edge  information  in  the  images  rather  than  directly  matching  pixels  in 
the  images  using  their  grey  levels.  Attempts  to  match  using  grey  level 
values  have  not  been  very  successful.  Crombie  [1975]  has  reported  anom¬ 
alous  matching  when  using  correlation  techniques  applied  directly  to  grey¬ 
scale  images.  Horn  and  Bachman  [1977]  used  a  hill  climbing  technique  on 
grey  scale  values,  but  their  method  requires  a  good  initial  approximate 
registration  and  is  also  computationally  expensive.  In  addition,  these 
problems  are  compounded  when  registration  of  images  from  different  sensors 
is  attempted.  For  example,  light-dark  reversal  is  possible  when  comparing 
a  panchromatic  image  with  an  infrared  image  of  the  same  scene.  Even  more 
of  a  problem,  are  radar  images  that  may  have  totally  black  shadows  where 
a  panchromatic  image  of  the  same  scene  shows  considerable  texture. 

One  argument  for  edge  matching  is  based  on  results  of  experiments  with 
human  beings.  People  seem  to  differentiate  regions  by  using  the  regions' 
boundaries.  John  Krauskopf  of  Bell  Laboratories  [Gilchrist  1979]  per¬ 
formed  experiments  on  human  beings  and  discovered  they  identify  the  color 
of  the  interior  of  a  region  using  contrast  information  at  the  boundary 
edges  of  the  region.  When  the  boundary  is  made  to  disappear,  by  artifi¬ 
cially  stabilizing  the  image,  the  perceived  color  of  the  region  changes 
to  the  color  of  the  surrounding  region. 

Another,  more  practical  argument,  is  that  edge-based  registration  of 
both  images  to  maps  and  of  images  from  different  sensors,  has  been  success¬ 
fully  performed.  Savol  [1978]  matched  radar  images  with  maps  using  edges. 
Hall  [1979]  matched  panchromatic  images  with  radar  images  using  correlation 
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techniques  applied  to  edges. 

This  section  deals  with  various  edge  based  feature  detectors. 

The  output  of  the  edge  detectors  is  a  binary  image  whose  pixels  with  value 

one,  are  the  corresponding  edge  points  in  the  original  image.  The  usefulness 
of  edges  for  registration  is  demonstrated  by  testing  the  results  of  some  of 
the  edge  detection  techniques.  This  was  accomplished  by  applying  the  edge 
detectors  to  two  dissimilar  images  of  the  same  scene  and  matching  the  resultant 
edge  images  using  the  simple  correlation  technique.  Although  the  application 
used  two  images,  it  just  as  easily  could  have  been  done  using  a  map. 
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3.2.1  First  Derivative  Methods 

One  simple  and  very  common  method  of  finding  edge  points  in  an  image 
is  to  apply  a  threshold  to  the  magnitude  of  the  gradient  (first  derivative) 
of  the  image.  There  are  many  ways  of  approximating  the  gradient  [Davis  1975] 
using  small  window  sizes.  Small  window  sized  methods  are  desirable 
because  they  are  computationally  fast;  some  methods  have  been  hardware  im¬ 
plemented.  Unfortunately,  they  also  are  vulnerable  to  noise  points.  More 
reliable  edge  points  can  be  found  by  using  larger  windows  [Rosenfeld  and 
Thurston,  1970;  Hayes  and  Rosenfeld,  1970]  but  the  ease  of  computation  de¬ 
creases  markedly.  In  all  cases  the  choice  of  the  appropriate  threshold  value 
has  to  be  determined. 

An  example  of  a  small-window  approximation  to  the  gradient  is  the  Sobel 
detector  which  is  defined  using  the  3x3  neighborhood  of  a  pixel  as  shown  in 
Figure  3-1. 

ABC 
D  E  F 
G  H  I 


Figure  3-1.  3x3  Neighborhood  of  Pixel  E. 

The  Sobel  magnitude  value  for  pixel  E  is  found  by 

IA+2B+C-G-2H-I  I+IA+2D+G-C-2F-I  I. 

The  Sobel  operator  was  applied  to  two  images  of  the  same  scene  but  taken 
from  different  sensors.  One  image  was  a  1024  x  1024  panchromatic  image  of  a 
scene  of  the  Canadian  border  (CBP52C).  The  other  image  was  the  corresponding 
infrared  image  (CBI52C) .  The  lower  left  512  x  512  quarter  of  the  panchromatic 
image  will  be  referred  to  as  image  BW-CB.  The  corresponding  quarter  of  the 
infrared  image  will  be  referred  to  as  IR-CB. 
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Using  a  threshold  of  166,  the  Sobel  operator  found  44,536  edge  points 
(17%)  in  image  BW-CB,  and  16%  in  image  IR-CB.  The  central  256  x  256  region 
of  detected  edge  points  from  the  images  were  correlated  and  part  of  the 
correlation  matrix  is  shown  in  Figure  3-2.  The  values  of  the  rest  of  the 
matrix  decrease  from  those  shown.  The  correlation  peak  occurs  at  the  dis¬ 
placement  (x  =  -4,  y  =  2). 
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Figure  3-2.  A  portion  of  the  correlation  matrix  of  the  central  regions 
of  Sobel  edge  points  from  images  BW-CB  and  IR-CB. 

The  experiment  was  repeated  using  a  threshold  value  of  260.  In  this 
case,  7%  of  BW-CB  pixels  and  6%  of  IR-CB  pixels  were  identified  as  edge  points. 
Unfortunately,  the  correlation  produced  five  peaks.  Two  high  peaks  occurred 
at  displacements  (3,  -5)  and  (-5,  -3),  while  there  were  three  other  peaks  at 
(-1,  1),  (0,  3),  and  (5,  3).  It  appears  that  setting  the  threshold  too  high 
will  cause  noise  points  to  be  preferentially  selected  making  correspondence 
unreliable.  Since  the  correct  transformation  is  known  to  be  (x=0,  v=0) .  the 
above  experiments  show  that  a  fairly  good  registration  transformation  can  be 
found  just  using  a  simple  edge  detector. 
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3.2.2  Marr  Edge  Detector 


Marr  [1979]  used  the  Laplacian  of  a  two-dimensional  Gaussian  distri¬ 
bution  convolved  with  an  image  as  an  edge  detector.  The  Laplacian  is  an 
approximation  to  the  second  derivative  of  the  image.  It  is  parameterized 
by  a  variable,  a,  so  that  the  Laplacian  can  be  made  smooth  and  band- 
limited  in  the  frequency  domain,  and  smooth  and  localized  in  the  spatial 
domain.  By  varying  a,  objects  within  a  determined  size  range  can  be 
detected  and  pinpointed  in  the  image. 

The  two  dimensional  Gaussian  is  defined  as 

1  -(x2+  y2) 

G(x’y)  =  2^  6  2o2 

and  its  Laplacian  is 

i  (x2+  y~) 

L(x,y)  =  G"  (x,y)  =  - — r-(x2+  y2  -  a2)e  2a2  . 

2ira° 

A  cross-section  of  L(x,y)  is  shown  in  figure  3-3 (a).  The  two- 
dimensional  L(x,y)  is  found  by  rotating  the  graphed  curve  about  the 
y-axis.  The  central  portion  of  L(x,y)  is  negative  and  is  2 o  wide. 

The  outer  portion  of  L(x,y)  is  positive.  The  peaks  of  L(x,y)  occur  at 
x  =  if3a  from  the  origin. 

A  cross-section  of  an  edge  is  shown  in  figure  3-3(b).  When  L(x,y) 
is  convolved  with  this  edge,  the  curve  in  figure  3-3 (c)  results.  Note 
that  where  this  curve  crosses  the  x-axis  is  where  the  edge  occurred  in 
figure  3—  3 (b ) .  Therefore  edges  are  found  by  locating  the  zero-crossings 
of  the  Laplacian  convolved  with  the  original  image. 

One  interesting  property  of  the  Marr  edge  detector  is  it  tends  to 
find  continuous  (gap-free)  edge  segments.  Most  other  edge  detection 
methods  do  not  do  this,  and  require  additional  processing  to  weed  out 


noise  points  and  fill  in  missing  pieces.  The  continuous  edge  segments 
can  then  be  easily  represented  as  chain  codes  [Freeman  1974]. 

The  Marr  edge  detector  with  c* 5.5  was  applied  to  the  BW-CB  and 
IR-CB  images.  The  outer  positive  portion  of  the  mask  was  chopped  off 
after  six  values  making  the  whole  mask  23  pixels  in  radius.  The  detector 
classified  17%  of  both  BW-CB  and  IR-CB  image  pixels  as  edge  points.  The 
central  256  x  256  portions  of  the  two  sets  of  edge  points  were  correlated 
and  part  of  the  resulting  correlation  matrix  is  shown  in  figure  3-4.  The 
correlation  peak  occurs  at  the  displacement  (x  =  -3,  y  =  3) .  Compared 
to  the  results  of  matching  using  the  Sobel  detector  (figure  3-2) ,  the 
Marr  detector  performed  much  better,  as  the  Marr  peak  value  has  twice  as 
many  matches  as  the  Sobel  correlation  peak  value. 

In  order  to  allow  for  corresponding  edges  in  two  images  to  be  off  by 
a  pixel,  it  may  be  desirable  to  thicken  the  detected  edges  [Hall  1979]. 

One  way  to  do  this  is  to  assign  a  weight  of  3  to  all  detected  edge  points, 
assign  a  weight  of  2  to  background  points  that  are  next  to  an  edge  point, 
and  a  weight  of  one  to  background  points  that  are  two  pixels  from  an  edge 
point.  All  remaining  background  points  are  assigned  a  weight  of  zero. 

The  thickened  edge  method  was  applied  to  the  detected  Marr  edge 
points  of  both  images,  and  the  central  256  x  256  portion  of  the  results 
were  again  correlated.  The  correlation  matrix,  in  this  case,  is  defined 
by 

3P4  334 

C(i,j)  =  l  E  T  (x,y)  *  T  (x  +  i,  y+j ) , 

X=129  y»129  1  *- 

where  T^(x,y)  and  T2(x,y)  are  the  values  of  the  thickened  edge  images  for 
BW-CB  and  IR-CB,  respectively. 
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Figure  3-5.  A  portion  of  the  correlation  matrix  of  the  central  regions 
of  thick  Marr  edge  points  from  images  BW-CB  and  IR-CB. 
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The  center  portion  of  the  correlation  matrix  is  shown  in  figure  3-5. 

The  peak  value  occurs  at  displacement  (x  =  -3,  y  =  3)  is  the  same  as  the 
result  obtained  from  the  unthickened  case.  However,  the  peak  value  in 
the  thickened  correlation  matrix  is  much  stronger  (since  170827/3 > 17089) . 
These  experiments  again  show  that  registration  using  edges  seems  to  work 
well. 

It  should  be  mentioned  that  the  Marr  detector  at  0=5.5  picked  up  tex¬ 
ture  elements  within  the  fields  and  woods  in  addition  to  the  boundaries 
between  woods  and  fields.  At  a  higher  value  for  a,  larger  regions  would 
be  picked  out.  Figure  3-6  shows  a  panchromatic  image  and  Figures  3-7 (a)- 
(c)  show  the  results  of  applying  the  Marr  detector  with  different  a  values. 
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Figure  3-6.  Image  used  to  demonstrate 
Marr  edge  detector. 
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3.2.3  Hough  Edge  Detector 

The  Hough  transformation  is  an  efficient  device  for  detecting  if  a  set  of 
high  contrast  points  are  organized  along  a  mathematical  curve  [Duda  and  Hart, 
1972].  The  simplest  mathematical  curve  and  the  most  important  one  for  detec¬ 
tion  of  man  made  structures  is  the  straight  line.  Only  two  parameters  are  re¬ 
quired  for  specification  of  a  given  line — in  polar  form  the  parameters  are  the 
direction  of  the  normal  to  the  line,  0,  and  the  distance  from  the  origin  to  the 
line,  r.  If  the  possible  line  directions  are  discretized  to  T  values  and  the 
possible  distances  from  the  origin  are  discretized  to  R  values  then  Hough  de¬ 
tection  is  logically  equivalent  to  a  matching  of  T-R  templates  to  the  high  gra¬ 
dient  points  [Stockman  and  Agrawala,  1977].  The  Hough  straight  line  detector 
is  an  excellent  device  for  detection  of  man  made  structures  in  imagery  as  proven 
by  various  L.N.K.  experiments  on  aerial  imagery. 

There  are  also  Hough  detectors  for  circles,  parabolas,  ellipses,  and  hyper¬ 
bolas.  LNK  has  performed  some  experiments  in  registration  using  the  Hough  cir¬ 
cle  detector,  [Stockman  and  Kopstein,  1979].  The  other  Hough  detectors  tend 
not  to  be  as  efficient  as  the  Hough  line  and  circle  detectors,  because  they 
require  estimation  of  too  many  parameters. 

Examples  of  the  result  of  applying  the  Hough  edge  detector  to  panchromatic 
images  are  shown  in  the  LNK  registration  procedure  section. 
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3.2.4  Relaxation 


It  is  possible  to  apply  the  technique  of  relaxation  to  edge  detection 
[Zucker  1977].  All  possible  edges  are  divided  into  groups  according  to 
their  orientations.  For  each  group,  the  probability  that  a  pixel  belongs 
to  an  edge  in  that  group  is  calculated,  plus  the  probability  that  the  pixel 
belonged  to  no  edge.  A  label  consisting  of  the  above  set  of  probabilities 
is  assigned  to  each  pixel.  Then,  through  a  parallel  iterative  process, 
the  labels  of  adjacent  pixels  are  compared  and  compatible  orientations 
cause  their  probabilities  to  strengthen  while  incompatible  orientations 
cause  their  probabilities  to  weaken.  The  process  stops  when  all  the  pro¬ 
babilities  in  all  the  labels  have  approached  their  limiting  values.  The 
initial  probabilities  can  be  found  by  using  the  simple  edge  detectors  dis¬ 
cussed  in  section  3.2.1. 

The  detected  edges  can  then  be  represented  by  a  chain  code.  The  re¬ 
laxation  procedure  tends  to  thicken  edges,  so  additional  processing  may 
be  necessary  to  compensate  for  this. 
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3.2.5.  Linking 

Most  of  the  edge  detection  methods  described  so  far  produce  only  edge 
points  (Sobel  detector)  or  edge  segments  (Hough  detector).  While  the  Marr 
detector  can  be  made  to  produce  relatively  long  edge  segments  it  is  computa¬ 
tionally  expensive.  It  is  possible  to  do  additional  processing  on  the  sim¬ 
pler  methods  to  try  to  produce  longer  edge  segments. 

One  such  method  takes  three  steps.  The  first  step  is  to  extract  all  high 
contrast  pixels.  A  gradient  operator  is  then  applied  to  each  high  contrast 
point  to  determine  the  gradient's  magnitude  and  direction  at  the  point.  In 
the  second  step,  a  small  neighborhood  about  each  high  contrast  point  is  spirally 
examined  to  find  the  best  (if  any)  continuing  high  contrast  point  in  the  forward 
and  backward  direction.  If  such  points  exist,  links  are  established  to  them. 
Finally  in  the  last  step,  the  edge  segments  are  extracted  from  the  chained  high 
contrast  points.  All  chains  less  than  a  certain  length  can  be  eliminated  as 
noise  edges.  Note,  the  linking  step  can  be  done  in  parallel  for  each  high  con¬ 
trast  point. 

LNK  has  applied  this  method  and  achieved  good  results  by  keeping  only  the 
top  5%  of  high  contrast  points.  Figure  3-8 (a)  shows  a  light  airplane  on  a 
darker  airfield.  The  result  of  applying  step  one  to  a  window  containing  the 
right  wing  tip,  i.e.  finding  the  high  contrast  points  and  applying  the  gradient 
operator  at  each  point  is  shown  in  Figure  3-8 (b).  The  result  of  applying  the 
second  step  of  linking  is  shown  in  Figure  3-8 (c).  Notice  that  some  of  the 
points  are  of  degree  three,  i.e.  they  are  at  the  junction  of  multiple  edge 
activity.  These  points  must  be  at  the  terminus  of  an  extracted  curve  segment 
because  they  cannot  relate  symmetrically  to  three  neighbors.  The  result  of 
applying  the  last  step  is  shown  in  Figure  3-8(d).  A  large  portion  of  the  wing 
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Figure  3-8(a). 
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Figure  3-8 (b), 


Gradient  direction  of  high  contrast 
points  of  right  airplane  wing. 
(Curve  detection  step  1.) 


Figure  3 -8(c).  Plot  of  all  forward  and  backward 
linking  relationships  among  high 
contrast  points  of  (b).  (Curve 
detection  step  2.) 
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boundary  was  successfully  extracted  along  with  two  edges  of  the  "USAF"  identi¬ 
fication  interior  to  the  wing  plus  two  edges  of  a  dark  streak  on  the  airfield 
below  the  plane. 

The  linking  of  edge  points  or  short  edge  segments  can  be  much  more  com¬ 
plicated.  Higher  level  processing  using  model-directed  techniques  can  be  used 
to  make  linking  decisions.  Geometric  or  topological  constraints  can  also  be 
used.  These  type  of  methods  produce  only  a  very  restricted  set  of  edges  and 
thus  are  not  very  general. 
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3 . 3  Point  Features 


Point  features  can  also  be  used  for  matching  and  in  some  ways,  may  be  more 
desirable  than  edge  features.  One  problem  with  edge  features,  is  that  the 
length  of  the  detected  edge  segments  is  likely  to  vary  among  images  of  the 
same  scene.  Since  it  is  difficult  to  tell  which  parts  of  the  edges  correspond, 
many  edge  matching  algorithms  use  the  end  points  of  the  edges.  But,  if  the 
endpoints  of  the  detected  edges  are  likely  to  vary,  this  makes  accurate  reg¬ 
istration  difficult.  Point  features,  on  the  other  hand,  do  not  have  this  prob¬ 
lem. 

Point  features  need  not  be  the  result  of  some  point  detectors.  For 
example,  the  intersection  of  two  lines  or  edges  determines  a  point.  Such 
point  features  are,  in  a  sense,  a  higher  order  feature  than  those  derived 
from  a  point  detector.  This  section  discusses  two  higher  order  point  fea¬ 
tures,  intersections  and  high  curvature  points,  and  two  point  detectors,  the 
Mo rave c  interest  operator  and  a  "special"  point  operator.  In  later  sections, 
methods  of  forming  edges  by  connecting  point  features  and  using  these  for 
registration,  are  discussed. 
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3.3.1  Intersections 


Intersections  between  two  or  more  lines  or  edges  can  provide  a  good  base 
for  registration.  The  topology  and  geometry  of  intersections  can  be  used  to 
classify  intersections  into  various  types.  For  example  the  number  of  lines  or 
edges  in  the  intersection,  the  angles  involved,  and  the  relative  grey  levels 
of  the  regions  between  the  lines  or  edges  can  all  be  used  to  type  an  inter¬ 
section.  The  typed  intersections  can  then  be  used  to  aid  matching. 

Since  most  edge  detectors  tend  to  be  unstable  at  intersections,  some  add¬ 
itional  processing  will  need  to  be  done  to  force  the  intersections  to  occur. 

It  is  even  possible  to  create  imaginary  intersections  as  the  surveyor  does. 

For  instance,  the  wall  of  a  building  can  be  extended  to  the  street  to  form  an 
intersection.  Matching  using  intersections  has  been  studied  [Zahn,  1974; 
Dudani,  1977;  Stockman  and  Kopstein,  1979]. 


3.3.2  High  Curvature  Points 


Not  all  images  will  contain  intersections.  Intersections  tend  to  stem 
from  man-made  structures  such  as  road  networks,  cultivated  fields,  or  build¬ 
ings.  In  images  with  no  man-made  features,  intersections  are  less  likely  to 
occur.  In  such  images,  high  curvature  points  on  curved  edges  can  be  used  as 
registration  features  instead. 

Once  long  continuous  edges  have  been  found,  they  easily  can  be  repre¬ 
sented  by  chain  codes  or  crack  codes.  Figure  3-9(a)  shows  a  small  region 
formed  by  '$'s  and  the  links  along  the  "cracks"  of  the  region.  Figure  3- 9(b) 
shows  how  these  links  can  be  represented  by  a  code.  The  crack  code  for  the 
links  of  the  region  in  Figure  3-9 (a)  is 
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Figure  3-9:  Crack  Code  Representation 


The  high  curvature  points  of  a  curve  represented  by  a  crack  code  car.  be 
defined  as  follows.  Let  the  1-curvature  of  a  point  P  on  a  curve  be  defined 
as  the  angle  between  the  line,  connecting  the  previous  point  on  the  curve  and 
P,  and  the  line,  connecting  P  and  the  next  point  on  the  curve.  The  K-curvature 
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can  then  be  defined  as  the  angle  formed  by  the  line,  connecting  the  Kcn  pre- 
vious  point  and  P,  and  the  line  connecting  P  and  the  succeeding  point. 

In  general,  the  curvature  measures  from  -179°  to  +179°  with  positive  values 
indicating  the  curve  is  bending  clockwise. 

For  example,  Figure  3-10  shows  a  portion  of  the  crack  code  of  a  curve. 

R  is  the  sixth  point  before  point  P  and  Q  is  the  sixth  point  after  point  P. 
The  K-curvature  for  point  P,  for  K=6,  is  the  angle  a  shown,  which  is  90°. 
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Figure  3-10.  An  example  of  K-curvature  for  point  P 
where  K  is  6.  The  K-curvature  is  a 
which  is  90° . 

A  technique  called  non-maximum  suppression  can  be  used  to  isolate  the 
high  curvature  points.  The  curvature  of  each  point  is  compared  with  its  K 
neighbors  on  both  sides.  If  the  curvature  of  the  point  is  more  positive  or 
more  negative  than  all  of  its  2K  neighbors,  then  it  is  a  high  curvature  point 
if  its  absolute  curvature  is  greater  than  90°, 

This  technique  was  applied  to  the  Marr  edges,  extracted  from  the  images 
BW-CB  and  IR-CB,  discussed  in  Section  3.2.2,  with  K=10.  Since  a  o=5.5  was 
used,  the  Marr  detector  picked  out  many  small  regions  (see  section  3.2.2)  so 
that  there  were  many  high  curvature  points.  A  higher  value  for  a  would  cause 


larger  regions  and  fewer  high  curvature  points;  thus  a  higher  value  for  K  would 
be  needed. 

The  high  curvature  points  of  the  two  images  were  correlated  and  part  of 
the  correlation  matrix  is  shown  in  Figure  3-11.  The  correlation  peak  occured 
at  displacement  (x=-2,  y=0)  which  is  closer  to  the  correct  transformation  than 
the  ones  obtained  just  using  the  edges.  However,  in  this  case,  the  full  512  x 
512  region  was  correlated  instead  of  the  central  256  x  256  region  as  before  and 
this  may  be  the  cause  of  the  difference. 
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Figure  3-11.  A  portion  of  the  correlation 

matrix  of  high  curvature  points 
from  the  images  BW-CB  and  IR-CB. 

The  cluster  in  figure  3-11  is  elongated  along  the  x-axis.  This  is  because 
there  are  several  ragged  horizontally-oriented  edges  in  the  images.  The  elonga¬ 
tion  could  then  have  been  caused  by  mismatching  a  high  curvature  point  with  a 
neighbor  of  the  corresponding  high  curvature  point  in  the  other  image.  All 
registration  techniques  have  problems  matching  regular  patterns  although  what 
is  considered  regular  depends  upon  the  technique  involved. 
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3.3.3  Ocher  point  features 


Moravec  [1977]  proposed  an  "Interest  Operator"  as  a  point  feature  detector. 
The  value  of  the  interest  operator  is  the  minimum  value  of  the  variance  of  the 
gray  s ',ale  values,  computed  in  four  directions  (horizontal,  vertical,  and  the 
two  diagonals)  over  a  small  window. 

The  imerest  operator  is  especially  sensitive  to  corners,  spots  and  un¬ 
fortunately,  noise.  The  operator,  of  a  particular  window  size,  would  be  applied 
to  the  whole  image  and  points  whose  interest  value  were  relatively  high  would 
be  selected.  Care  in  picking  the  threshold  would  be  needed  so  that  noisy  points 
are  not  the  only  points  chosen. 

It  is  possible  to  define  a  set  of  "special"  points  for  each  image  by  finding 
the  "special"  point  in  each  row,  column,  and  diagonal  in  the  image.  These 
special  points  are  found  by  applying  any  long  one-dimensional  detector  to  each 
row,  column,  and  diagonalj  and  selecting  the  point  with  the  highest  detector 
response. 

There  are  several  advantages  to  this  method.  The  first  is  that,  only  a  few 
points  will  be  selected  which  cuts  down  the  computation  of  matching.  Second,  the 
points  are  automatically  distributed  at  least  somewhat  throughout  the  image,  since 
each  row,  column,  and  diagonal  will  have  a  special  point.  This  method  would  work 
best,  either  matching  images  from  the  same  sensor,  or  matching  images  to  maps. 
Note,  the  map  need  only  store  a  small  number  of  points. 

The  one  dimensional  detector  could  be  an  edge  detector,  line  detector,  tex¬ 
ture  detector,  and  so  on.  An  experiment  using  a  one-dimensional  edge  detector 
was  performed  on  images  BW-CB  and  IR-CB,  described  in  Section  3.2.1.  For  each 
row  and  column,  (the  diagonals  were  not  used),  the  point  which  had  the  greatest 
difference  between  the  sum  of  the  preceding  50  points  and  the  sum  of  the  fol- 
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lowing  50  points,  was  chosen  as  the  special  point.  The  special  points  were 
correlated  and  a  portion  of  the  correlation  matrix  is  shown  in  Figure  3-12. 
The  cluster  is  again  elongated  along  the  horizontal  as  in  Figure  3-11  and 
probably  for  the  same  reason,  i.e.  mismatching  of  neighboring  special  points 
along  a  prominent  edge.  The  peak  occurs  at  (x=-3,  y=0)  which  is  different 
from  those  obtained  before. 
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3.4  Region  Features 
3.4.1  Invariant  Moments 


Features  on  images  or  parts  of  images  which  are  invariant  under  rotation, 
translation,  and  scale  changes  can  be  used  for  scene  matching.  It  is  possible 
to  define  such  invariant  features  using  the  moments  of  the  image  gray  scale 
function.  The  problems  with  moments  are  (1)  digitization  may  cause  the 
moments  to  no  longer  be  invariant  and  (2)  finding  an  easily  computable  set 
of  moments  that  can  adequately  characterize  an  image. 

Wong  and  Hall  [1978]  have  devised  a  set  of  seven  invariant  moments  which 
they  claim  are  adequate  for  matching  two  images.  Let  the  gray  scale  function 
be  denoted  f(x,y),  then  the  pq1^  central  moment  of  the  image  is  defined  by 

np9  =  x  y  (x-^p(y-y)qf  <x«y). 

where  x  and  y  are  the  means  of  x  and  y,  respectively,  over  the  whole  image. 

The  normalized  central  moments  pq  can  then  be  defined  by 
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where  y  =  (p  +  q)/2. 

The  set  of  seven  moments  which  are  invariant  under  translation,  rotation, 
and  scale  are  then  defined  by 

M1  =  B20  +  B02’ 

M2  '  <B20  '  B02>2  +  4BU2- 

M3  '  <B30  '  3B12)2  +  <3B21  +  b03>2- 

M4  *  «30  +  B12>2  +  <B21  +  V2’ 

M5  ’  <B30  +3B12)(B30  +  B12)l(B30  +  B12)2-3(B21  +  B03>21 


+  ( 3iJ21  b03^B21  +  B03^3^B30  +  B12^  "  ^B21  +  B03^  ^ 


/ 
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M6  =  (B20  -  B02)[(B30  +  B12)2  “  (B21  +  B03)2]  +  4B11(B30  +  B12> 

( B2i  Bq3)  * 

M?  =  (3B12  -  3B30)(B30  +  B12)[(B30  +  B12)2  -  3 (B21  +  BQ3) 2 ] 

+  (3B21  -  BQ3)x(B21  +  B03)[3(B30  +  B12)2  -  (B21  +  B03)2]. 

There  are  many  possible  ways  of  using  the  invariant  moments.  The  original 
image  could  be  subdivided  into  smaller  subimages  or  regions  and  each  subimage 
or  region  could  be  represented  by  its  own  set  of  invariant  moments,  and  these 
could  be  used  to  correlate  with  other  images  or  a  map.  Section  3.4.2  presents 
some  methods  of  region  segmentation  and  section  4.1.3  discusses  a  method  of 
finding  sub images. 
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3.4.2  Region  Segmentation 


A  segmentation  of  an  image  is  a  partition  of  the  image  into  disjoint 
subsets  whose  union  is  the  entire  image.  For  surveys  of  segmentation  see 
[Riseman  and  Arbib  1977,  Kanade  1980,  Zucker  1976,  Pavlidis  1977].  Many 
segmentation  procedures  produce  an  initial  segmentation  and  then  apply  an 
iterative  procedure  such  as  merging  or  splitting  to  obtain  an  improved  seg¬ 
mentation.  Two  basic  approaches  to  finding  an  initial  segmentation  are  first 
locating  boundaries  or  first  locating  pixels  with  similar  feature  values. 

Boundary  detection  generally  consists  of  edge  detection  followed  by 
linking  of  edges  into  closed  boundaries  as  discussed  in  section  3.2.5.  Edge 
detection  procedures  which  form  closed  edges  such  as  the  Marr  detector  (sec. 
3.2.2)  and  relaxation  labelling  (section  3.2.4)  provide  a  segmentation  direct¬ 
ly.  Edge  detection  in  textured  images  is  a  difficult  task  which  depends  heav¬ 
ily  on  the  relative  scale  of  the  texture  and  the  regions'  sizes. 

Construction  of  regions  based  on  similarity  of  pixels  or  pixel  neighbor¬ 
hoods  requires  feature  extraction.  Features  commonly  measured  include  average 
gray  level  in  a  neighborhood,  variance,  average  edge  content  per  unit  area, 
average  orientation  of  local  edges  and  average  spot  size  of  uniform  contiguous 
areas.  Texture  measures,  such  as  co-occurrence  matrices  and  Fourier  transform 
ring  data  can  be  used  as  feature  measures  for  larger  areas.  Threshold  tech¬ 
niques  [Pavlidis  1977,  Price  1976]  can  be  used  with  these  features  to  provide 
an  initial  segmentation. 

Milgram  [1978]  describes  a  procedure  for  region  construction  using  evi¬ 
dence  from  several  sources  such  as  edge  information  and  pixel  feature  values. 
The  algorithm  selects  the  contours  at  different  thresholds  according  to  the 
support  of  the  edge  data  along  the  contours.  Zucker  [  1979  ]  gives  a  relaxa¬ 

tion  technique  for  constructing  regions  from  primitive  edges.  This  scheme 
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allows  points  to  be  considered  as  interior  or  boundary  points  of  a  region. 

The  relaxation  process  allows  edge  segments  separating  regions  to  prosper 
as  region  points  and  edge  points  reinforce  themselves. 

Pavlidis  describes  a  general  class  of  split-merge  algorithms  using  the 
notion  of  region  adjacency.  Given  a  criterion  for  deciding  if  a  single  region 
should  be  split  and  a  criterion  for  deciding  whether  two  adjacent  regions 
should  be  merged,  the  procedure  is  as  follows: 

1)  Split  each  region  which  should  be  split  according  to  the  splitting 
criterion.  Continue  this  until  no  further  regions  satisfy  the  split¬ 
ting  criterion. 

2)  If  step  3  has  not  been  executed  yet  then  keep  going,  else  if  no 
merges  occurred  in  the  last  execution  of  step  3,  then  stop. 

3)  Merge  any  two  adjacent  regions  which  should  be  merged  according  to 
the  merging  criterion.  Continue  this  step  until  no  adjacent  regions 
satisfy  the  merging  criterion. 
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4.  Registration 

The  goal  of  registration  is  to  find  a  correspondence  of  pixels  between 
two  images  or  between  an  image  and  a  map.  In  this  section  we  survey  some 
existing  techniques  for  matching.  Before  proceeding  with  the  description  of 
current  methods  we  will  attempt  to  give  a  general  framework  for  discussing 
matching  problems. 

Matching  involves  a  comparison  of  representations.  We  may  classify  match¬ 
ing  algorithms  by  the  type  of  representation  we  are  comparing.  The  representa¬ 
tion  may  be  the  original  image,  a  hierarchy  of  images  of  varying  sizes  and  re¬ 
solutions,  an  image  derived  from  the  original  image,  a  graph,  or  a  list  of 
points,  edges  or  regions  together  with  a  set  of  descriptors.  Derived  images  may 
consist  of  regions  of  constant  intensity,  region  boundaries,  edges  or  isolated 
points.  Matching  techniques  depend  upon  the  representation  selected.  Correla¬ 
tion  techniques  can  be  applied  to  images,  graph  isomorphism  techniques  may  be 
applied  to  graphs,  relaxation  procedures  may  be  applied  to  lists,  and  clustering 
may  be  applied  to  points,  edges,  and  abstract  vectors  and  triangles.  A  hierarchy 
of  derived  images  may  be  used  to  locate  subregions  of  an  image  which  are  likely 
to  provide  a  quick  accurate  match. 

i 

The  representations  and  matching  procedures  we  have  referred  to  admit  nu¬ 
merous  modifications  and  combinations.  Since  no  firm  guidelines  exist  for 
selecting  matching  procedures,  we  indicate  some  relevant  considerations.  For 
all  methods  other  than  the  standard  correlation  of  gray-scale  imagery,  the  in¬ 
expensive  extraction  of  good  features  is  a  primary  concern.  The  cost  of  fea¬ 
ture  extraction  must  be  balanced  against  the  potential  use  of  the  features. 

Region  features,  while  often  costly  to  obtain,  can  be  vital  in  later  phases  of 
image  interpretation.  On  the  other  hand,  isolated  interest  points,  such  as 
points  of  high  curvature,  can  be  inexpensive  to  obtain  but  serve  no  purpose 
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other  than  to  determine  registration.  A  second  major  consideration  is  the 
extent  to  which  knowledge  of  the  scene  contents,  as  might  be  obtained  from  a 
geographic  data  base,  should  be  employed  in  the  matching  procedure. 


I 

I 

I 

I 

I 
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4  1  Correlation-Based  Techniques 

Registration  can  be  done  using  various  correlation-based  techniques. 

These  techniques  can  be  applied  to  the  original  gray-scale  imagery  or  to  de¬ 
rived  pixel  images  such  as  images  which  contain  point  features  or  edges,  or 
can  be  applied  to  lists  of  derived  feature  vectors. 

Because  of  the  high  cost  and  ineffectiveness  of  basic  correlation,  other 
methods  have  been  derived.  This  section  will  present  the  basic  correlation 
method,  a  sequential  correlation  method  and  a  hierarchical  correlation  method. 

While  correlation  can  be  used  to  determine  the  approximate  global  trans¬ 
lation  for  registration,  the  cost  of  determining  the  proper  rotation  is  general¬ 
ly  prohibitive  and  no  provisions  are  made  for  accommodating  scale  change  or 
perspective.  Correlation  techniques  should  work  with  images  and  maps,  although 
little  work  has  been  done  on  this.  Edge-based  correlation,  whether  real  images 
or  abstract,  give  the  best  chance  for  registration. 
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4.1.1  Basic  Correlation 


The  correlation  function  R(u,v),  for  a  pixel  image,  at  a  point  (u,v) 
of  a  JxK  window  W  with  an  MxN  image  is  defined  as 


R(u,v) 


£ 


i=l 


K  J 

£  £  S (i, j )W(i-u, j-v) 

j=l  i-1 _ 


s2(i.j)i1'2[£  £ 

j=l  i=l 


W2(i-u,j-v)]1/2 


2 

since  W  (i-u,j-v)  is  constant  for  a  fixed  W,  it  suffices  to  compute 

S(i, j )W(i-u, j-v 


R(u,v) 


i=l 


-v) 


S  S2  (i,  j  )  l1^2 
J-l  i=l 


The  point  (u,v),  in  the  search  image,  at  which  R(u,v)  is  maximized  gives  the 
optimal  location  of  the  window.  Thus  our  optimal  match  is  given  by  the 
(u0,v0)  such  that 

R(u0,v0)  >  R(u, v)  for  1  <  u  <  (M-J+l) 

1  <  v  <  (N-K+l) 

The  above  pixel  based  matching  procedure  can  be  applied  to  original  gray 
scale  images,  edge  images  and  point  images.  These  methods  can  be  easily 
adapted  to  the  case  where  the  image  is  described  by  a  feature  vector,  such  as 
the  vector  of  invariant  moments  of  an  image  (see  section  3.4.1).  If  we  de¬ 
note  the  seven  invariant  moments  of  image  1  by  M^,  ...,  and  the  invariant 
moments  of  image  2  by  N^,  ...»  Ny  then  we  define  their  correlation  by: 


Invariant  moment  matching  is  in  theory  invariant  under  translation,  rotation, 
and  scale  change,  though  in  practice  good  results  of  matching  radar  to  optical 
images  have  only  been  obtained  for  angle  of  up  to  45°  and  scale  changes  of  a 
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factor  of  2,  [Wong  and  Hall,  1978].  Invariant  moment  matching  may  be  used 
as  part  of  a  hierarchical  procedure  discussed  later. 

Correlation-based  techniques  can  be  applied  to  lists  of  features.  Of 
particular  importance  is  the  case  where  a  feature  such  as  a  straight  line  or 
a  triangle  can  be  assigned  to  certain  points  in  the  image.  Given  such  an 
assignment  for  two  images  and  a  transformation  from  one  image  to  the  other, 
we  would  like  to  use  the  features  to  measure  the  quality  of  the  transforma¬ 
tion  as  a  registration.  The  measure  will,  of  course,  depend  upon  the  nature 
of  the  feature.  In  the  comparison  of  straight  lines  from  two  images,  we  can 
define  two  lines,  one  from  each  image,  to  be  similar  if  the  sums  of  the  squares 
of  the  distances  of  their  endpoints  are  close  together.  Since  a  triangle  can 
be  thought  of  as  three  straight  lines,  we  can  extend  the  previous  definition 
to  the  case  of  triangles  in  a  straight  forward  manner. 

The  correlation  function  of  two  line  images  can  be  defined  as  follows: 

Let  a  be  a  vector  giving  the  parameters  of  a  transformation  mapping  one 
image  into  another,  e.g.  translation,  rotation  and  scale  parameters. 

The  correlation  function  assigns,  to  each  a,  a  number  f(a),  measuring 
the  quality  of  the  transformation  for  registration. 

To  each  line,  a,  in  image  one  which  has  a  line,  b,  in  image  two  whose  endpoints 
are  each  within  a  fixed  value,  e,  of  the  endpoints  of  a,  we  assign  the  value 
Da  which  is  the  sum  of  the  distances  between  the  endpoints  of  the  two  lines. 

Let  A={a^,  ....  a^}  be  the  set  of  all  edges  in  image  one  for  which  edges  b,  as 
above,  exist.  Define  the  correlation  f (a)  by 

f(a)  =  “g-  (  Yi  D*  +  D  if  WO 

a.e  i 
=  h  otherwise 

where  h  is  a  large  positive  constant  and  g  is  an  integer  parameter  to  be  set. 
The  a  which  minimizes  f  is  said  to  give  the  best  correlation. 
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4. 1.1.1  Correlating  Abstract  Edge  Lists 

An  example  of  correlating  feature  vectors  where  the  features  are  edges 
is  presented.  Suppose  edges  of  Marr  high  curvature  points,  discussed  in  Sec¬ 
tion  3.3.2,  are  formed  by  connecting  successive  high  curvature  points.  The 
resulting  edge  list  does  not  contain  real  edges  but,  what  LNK  calls,  abstract 
edges.  This  process  was  applied  to  the  images  BW-CB  and  IR-CB  and  the  result¬ 
ing  abstract  edge  lists  were  correlated.  The  edges  were  said  to  correspond 
when  the  orientation  of  the  edges  were  within  ±30°  and  the  lengths  of  the  edges 
were  within  one  pixel.  Part  of  the  correlation  matrix  is  shown  in  Figure  4-1. 
The  peak  occurred  at  displacement  (x=-2,  y=0)  which  is  better  than  the  straight 
correlation  of  the  high  curvature  points.  Note,  the  peak  is  much  better  formed 
than  the  peaks  in  the  pixel  based  matching  correlations. 


45 


Figure  4-1.  Part  of  the  correlation  matrix  of  abstract  edge 
of  high  curvature  points. 


4. 1.1. 2  Correlating  Abstract  Triangle  Lists 

The  example  of  the  previous  section  was  repeated  but  using  abstract  tri¬ 
angles  formed  by  connecting  three  successive  high  curvature  points.  Figure 
4-2  shows  part  of  the  resulting  correlation  matrix.  The  triangles  were  said 
to  correspond  if  their  size  and  orientation  were  similar.  Again,  the  peak 
occurred  at  (x=-2,  y=0) .  These  two  experiments  show  that  correlation  using 
abstract  features  can  work  much  better  than  the  basic  pixel-based  correlation. 


4.1.2  Sequential  Matching  Methods 


One  method  for  reducing  the  amount  of  computation  in  performing  matching  is 

to  do  the  matching  sequentially  [Barnea  &  Silverman  1972].  Suppose  there  is  an 

LxL  array  of  pixels,  W,  called  a  window  and  an  MxM  array  of  pixels,  S,  called  a 

search  area  with  M<L.  The  aim  of  the  procedure  is  to  find  that  MxM  subarray  of 

S  which  correlates  best  with  W.  Let  denote  the  MxM  subimage  of  S  with  the 

reference  point  (i,j)  at  the  upper  left  comer.  Define  coordinates  in  the  array 

S1, j  by 
M  7 


S^’J  (n,m)  =  S (i  +n-l,  j+m-1) 
for  1  <  n,  m  <  M 


1  <  i,  j  <  L  -  M  +  1. 

The  window  W  is  compared  with  its  subwindow  SjJj’J  by  comparing  a  subset  of 

its  corresponding  pixel  pairs,  defined  by  <S^j,J(n,m),  W(n,m)>.  The  possible 

2  2  ? 

number  of  comparisons  required  is  M  (L-M+l)  ,  since  there  are  (L-M+l)  windows 

o  2 

each  of  size  M*-.  The  sequential  matching  method  requires  examining  all  (L-M+l) 

2 

windows  but  not,  in  general,  all  of  the  M  pixel  pairs.  A  measure  of  dissimilar¬ 
ity  between  pixel  pairs  is  defined  so  that  matching  of  windows  can  be  accomplished 
by  accumulating  the  dissimilarity  of  pixel  pairs  until  a  threshold  is  reached. 

In  order  to  do  this,  the  pixel  pairs  must  be  ordered  in  some  way. 

2 

One  procedure  for  ordering  the  M  pixel  pairs  is  by  randomly  generating  a 

2 

permutation  of  the  numbers  1,  ...,  M  and  examining  the  pairs  in  the  order  speci¬ 
fied  by  the  permutation.  Then  two  measures  of  dissimilarity,  or  error,  between  the 
components  of  the  pixel  pair  S^’J(p,q)  and  W(p,q)  can  be  defined.  They  are  e’, 
the  unnormalized  measure,  and  c,  the  normalized  measure,  given  by 
s'(i,j,p,q)  5  |S*,:i(p,q)  -  W(p,q)j, 
e(i,j>p,q)  =  |s^(p,q)  -  S (i ,  j  )  -  W(p,q)  +  w]  , 

M 
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where 


M  M 


^  =  rj2~  2  W(n,m), 

M  n=l  m=l 


MM 


=  M2  S  SM 


n=l  m=l 

The  errors  e  or  e '  are  accumulated  by  comparing  the  pixel  pairs  in  the  order 
specified  by  the  permutation,  until  a  given  threshold  is  reached.  The  number  of 
pixel  pairs,  of  S^’-*  and  W,  compared  before  the  threshold  is  surpassed  is  stored 
in  an  array  F  at  the  location  F(i,j).  The  points  (i,j),  for  which  F(i,j)  is 
large  correspond  to  potential  window  matches.  The  potential  window  matches  can 
then  be  used  by  other,  more  exhaustive,  matching  procedures  as  initial  guesses. 

Variations  of  the  sequential  method  have  been  described  [Barnea  &  Silverman, 
1972;  Ramapriyan,  1976].  For  example,  the  threshold  may  be  given  as  a  function 
of  the  rate  of  change  of  accumulated  error  with  respect  to  number  of  point  pairs 
examined  so  far.  In  addition,  other  methods  for  ordering  the  pairs  are  possible. 

As  with  basic  correlation,  sequential  correlation  is  designed  to  determine 
a  translation  for  registration.  Rotational  correction  is  still  expensive  and 
scale  changes  are  not  considered.  Hall  describes  sequential  correlation  of 
edge  images  which  could  also  be  applied  to  the  registration  of  an  image  to  a 
map.  This  edge  image  correlation  may  be  applied  to  images  from  dissimilar  sen- 
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4.1.3  Hierarchical  Scene  Matching 


Hierarchical  scene  matching  [Wong  &  Hall  1978]  is  a  technique  for  register¬ 
ing  one  image  with  another,  using  a  sequence  of  images  derived  from  the  original 
images,  but  varying  in  size  and  resolution.  Each  successive  set  of  images  are 
smaller  in  size  and  have  lower  resolution.  Let  level  k=0  correspond  to  the 
original  image  which  is  assumed  to  be  square.  The  (k+1)—  set  of  images  is 
found,  from  the  k—  set,  by  partitioning  each  k-level  image  into  four  equal  size 
square  windows,  applying  a  low-pass  filter  to  each  window  and  then  sampling  at 
1/2  the  sample  rate  of  the  previous  level.  This  gives  a  set  of  smaller  windows 
at  a  lower  resolution. 

At  some  point,  this  process  is  stopped  and  the  highest  level  sub images  of 
the  two  original  images  can  be  compared.  Many  registration  methods  could  be 
used  for  the  comparisons,  such  as  the  previously  discussed  sequential  method. 

The  hierarchical  technique  need  not  be  applied  to  the  original  gray-scale 
image.  Instead  binary  edge  images  could  be  extracted  and  the  hierarchical 
process  could  be  applied  to  them,  for  example.  Or,  each  sub image  could  be 
represented  by  a  set  of  the  seven  invariant  moments  presented  in  section  3.4.1 
and  these  could  be  correlated. 

As  in  Section  4.1.2,  hierarchical  correlation  is  designed  to  determine  a 
registration  translation.  Due  to  the  computational  savings  offered  by  this 
method  over  ordinary  correlation,  it  may  be  possible  to  take  into  account  rota¬ 
tion  and  scale  changes.  Hall  [1979]  describes  an  application  of  this  method  to 
matching  optical  to  radar  images.  In  this  work  he  performs  the  hierarchical 
matching  on  edge  images  formed  from  the  original  images.  This  method  could  be 
adapted  to  the  matching  of  an  image  to  a  map. 


4.2.  Basic  L.N.K.  Registration  Technique 


This  section  presents  an  overview  of  the  basic  LNK  registration  procedure. 
A  more  detailed  description  can  be  found  in  Appendix  A.  The  goal  of  registra¬ 
tion  is  to  use  automatically  extracted  features  from  two  images  (or  an  image 
and  a  map)  to  find  the  best  global  transformation  that  maps  one  image  to  the 
other. 

The  LNK  procedure  is  able  to  achieve  this  even  when  many  local  mismatches 
of  features  occur.  The  three  basic  steps  of  the  method  are 

1)  Primitive  features,  such  as  line  segments,  edge  segments,  inter¬ 
sections,  or  high  curvature  points  are  automatically  extracted. 

These  features  may  be  parameterized,  such  as  by  position,  by 
length,  or  by  orientation. 

2)  Assume  all  features  of  one  type  can  correspond  to  one  another. 

For  example,  an  intersection  of  three  lines  in  one  image  can 
correspond  to  any  three-line  intersection  in  the  second  image. 

However,  only  one  of  these  correspondences  is  the  true  one. 

Now,  for  each  of  the  possible  correspondences,  find  the  trans¬ 
lation  and  rotation  that  maps  one  feature  to  the  other.  Let 
that  translation  and  rotation  transformation  be  represented  by 
the  triple  (A0,  Ax,  Ay).  Place  a  unit  of  weight  in  the  bin  in 
the  three  dimensional  histogram  that  represents  (AO,  Ax,  Ay). 

This  process  is  repeated  until  all  possible  feature  correspon¬ 
dences  have  been  transformed. 

3)  Locate  any  prominent  clusters  of  (A0,  Ax,  Ay)  in  the  histogram. 

Each  cluster  represents  a  set  of  features  in  one  image  that 
could  be  mapped  to  corresponding  features  in  the  other  image  by 
one  particular  (AO,  Ax,  Ay).  The  best  global  transformation  is 
defined  by  the  (A3,  Ax,  Ay)  of  the  largest  cluster.  This  trans¬ 
formation  provides  the  largest  number  of  local  correspondences. 

This  method  works  because  the  correct  transformation  will  show  up  as  a 
large  cluster,  while  the  wrong  transformations  will  tend  to  be  distributed 
randomly  throughout  the  histogram.  The  procedure  may  also  be  performed  it¬ 
eratively  by  making  the  bin  sizes  smaller  and  smaller. 

As  an  example,  suppose  an  image  can  be  represented  by  the  four  directed 
edge  segments  shown  in  Figure  4-3 (a).  The  direction  of  the  arrow  is  deter¬ 
mined  by  making  the  region  on  the  right  of  the  edge  segment  the  darker  one. 
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Figure  4-3.  Example  of  basic  LNK  registration 
technique.  Image  edge  elements  in 
(a)  need  to  be  rotated  45°  and  then 
translated  (4.5,  -2.0)  to  be  trans¬ 
formed  into  corresponding  map  edge 
elements  in  (b) . 


Suppose  the  map  or  second  image  contains  the  edge  segments  in  Figure  4-3 (b). 
There  are  16  possible  ways  that  the  edge  segments  in  (a)  can  be  paired  with 
the  edge  segments  in  (b).  Four  of  the  16  pairings  yield  a  consistent  inter¬ 
pretation — rotate  by  0=4.5°  and  translate  by  (4.5,  2.0).  The  triple  (0*45°, 
Ax=4.5,  Ay=2.0)  forms  a  cluster  in  the  three-dimensional  histogram  while  the 
triples  from  the  incorrect  pairings  are  sparsely  distributed.  Table  4-1  shows 
the  16  possible  transformations.  The  correct  transformations  are  indicated  by 
an  In  actual  cases,  there  will  be  many  more  than  4  primitive  features 

and  not  all  pairings  will  be  possible  (i.e.  due  to  size,  shape,  or  type  differ¬ 
ences)  so  the  cluster  of  correct  transformations  should  be  even  more  prominent. 

As  stated  in  step  1  of  the  procedure,  many  features  could  be  used  to  per¬ 
form  the  transformation.  Not  all  features  are  equally  desirable.  For  example, 
the  lengths  of  edge  segments  or  line  segments  usually  can  not  be  accurately 
determined.  Small  perspective  changes  can  alter  the  curvature  of  high  curvature 
points  considerably.  LNK  feels  that  intersections  would  produce  more  accurate 
transformations  as  they  can  be  determined  more  accurately.  LNK  has  performed 
enough  experiments  to  show  that  even  if  up  to  90 Z  of  the  detected  features  in 
the  image  do  not  match  with  any  features  in  the  map,  the  correct  transformation 
will  still  form  enough  of  a  cluster  to  be  detected. 
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4.2.1  Example 


The  basic  LNK  registration  technique  was  applied  to  the  image  in  Fig¬ 
ure  4-4.  The  straight  line  Hough  detector  was  applied  to  this  image  and  the 
result  is  shown  in  Figure  4-5.  Another  view  of  the  naval  base,  shown  in  Figure 
4-6,  was  used  to  create  a  map  of  edge  segments,  shown  in  Figure  4-7.  Three 
ground  control  points  were  used  in  each  image  to  establish  the  approximate 
registration  transformation,  (0=328°,  Ax=-141,  Ay=9) . 

A  hierarchical  clustering  technique  was  used  to  determine  the  results.  In 
this  method,  clustering  is  first  done  on  0  and  then  in  (Ax,  Ay)  space  given  a 
fixed  0.  The  result  is  shown  in  Table  4-2.  The  strongest  transformation, 
(0=330°,  Ax=-142,  Ay=8)  aligned  19  out  of  200  image  edge  segments  with  12  out 
of  22  map  edge  segments.  The  strongest  contender  (0=237°,  Ax=453,  Ay=-433) 
aligned  only  6  image  edge  segments  with  2  map  edge  segments.  The  results  would 
probably  be  enhanced  if  a  more  comprehensive  map  had  been  constructed. 


Table  4-2 


Summary  of  Results  Using  the 
Basic  LNK  Registration  Technique 


1  0  cluster 

(Ax, 

Ay)  cluster 

index 

center /radius /strength 

index 

center /radius /strength 

cluster 

1 

330°/3°/176 

cluster  11 

(-71,  -215) /10/  2.35 

cluster  12 

(-142,  8)  /10/12.85 

cluster  13 

(195,  135)  /10/  5.04 

cluster 

2 

237°/3°/152 

cluster  21 

(453,  -433)/20/  5.14 

cluster  22 

(345,  78)  /10/  3.48 

cluster  23 

(276,  146)  /10/  3.86 

cluster 

3 

61°/3°/139 

cluster  31 
* 

* 

(-22,  -406)/20 /  1.28 

*  indicates  no  viable  alternative  cluster 
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Figure  4-5.  Hough  detections  from  UNB1  Figure  4-4  image 
with  2°  directional  resolution. 
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4.2.2  Full  RST  Transformation 

The  position  and  orientation  of  primitive  features  such  as  intersections, 
high  curvature  points,  or  lines  can  accurately  be  determined,  but  it  is  diffi¬ 
cult  to  determine  their  sizes,  such  as  length.  The  basic  LNK  registration 
technique  uses  the  position  and  orientation  of  the  primitive  features  to  find 
the  rotation  and  translation  necessary  to  register  two  images  or  to  register 
an  image  to  a  map. 

If  the  size  could  be  included,  then  it  would  be  possible  to  calculate 
the  four  parameter  transformation  of  rotation,  translation  and  scaling.  This 
section  presents  a  method  for  extending  the  procedure  to  account  for  scaling. 

In  order  to  accomplish  this,  abstract  vectors  or  edges  whose  size  can  accurately 
be  determined  are  introduced. 

To  achieve  scaling,  abstract  vectors  or  edges  can  be  formed  by  spanning 
pairs  of  point  structures.  For  example,  abstract  edges  can  be  formed  by  con¬ 
necting  pairs  of  high  curvature  points.  Abstract  vectors  could  be  formed  by 
using  an  intersection  point  as  the  vector  tail  and  a  high  curvature  point  as 
the  vector  head.  There  are  many  ways  of  forming  the  abstract  edges  or  vectors. 

The  registration  procedure  is  similar  to  the  basic  one.  Instead  of  the 
triples  (0,  lx.  Ay),  there  are  four  parameters  (0,  Ax,  Ay,  As),  where  As  is 
the  scaling  parameter. 

The  three  registration  steps  are  now:  (Appendix  A  provides  a  detailed 
description) 

1)  Primitive  point  features  (intersections,  high  curvature  points, 
etc)  are  automatically  extracted.  The  abstract  vectors  are 
created  by  pairing  the  primitive  point  features.  Not  all  pos¬ 
sible  pairs  need  be  formed  as  that  would  result  in  a  great  deal 
of  computation. 

2)  Assume  all  features  of  one  type  can  correspond  to  one  another. 

That  is,  a  vector  from  an  intersection  of  three  lines  to  an 
intersection  of  four  lines  in  the  first  image  can  correspond 
to  any  3-intersection  to  4- intersection  vector  extracted  from 
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the  map  (or  second  image).  For  each  possible  correspondence,  find 
the  4-parameter  transformation  (Ac,  Ax,  Ay,  As)  that  maps  one  vector 
to  the  other.  Place  a  unit  of  weight  in  the  bin  in  the  four  dimen¬ 
sional  histogram  that  represents  the  (AO,  Ax,  Ay,  As)  found. 

3)  Locate  the  best  cluster  in  the  histogram.  The  (AO,  Ax,  Ay,  As)  of 
that  cluster  is  the  best  global  transformation  as  it  provided  the 
largest  number  of  local  correspondences. 

As  an  example  of  this  method,  suppose  an  image  and  map  are  represented  in 
terms  of  vectors  v.  connecting  intersection  points.  The  intersection  points  are 
of  four  types;  'L' ,  'T',  'X',  or  'Y'.  The  rules  for  pairing  points  to  form 
vectors  may  be  arbitrary;  for  example,  'T'  points  are  joined  with  'L'  points. 

The  purpose  of  such  rules  is  to  control  the  combinatorics.  In  Figure  4-8,  five 
vectors  represent  the  map  and  four  vectors  represent  the  image. 

Given  any  map  vector  Vj  all  possible  matching  image  vectors  are  con¬ 
sidered.  Each  possible  pairing  (v.,v.)  results  in  an  RS&T  transformation 

1  ^ 

mapping  image  vector  onto  map  vector  v..  The  transformation  is  specified 
by  a  quad  of  parameters  (A3,  Ax,  Ay,  As)  where  AO  is  the  angle  of  rotation, 

As  is  the  scale  factor,  and  Ax  and  Ay  are  the  x  and  y  translations  respectively. 
The  pair  is  discarded  without  producing  a  quad  if  the  tips  or  tails  of 

the  vectors  and  disagree  in  type.  In  the  example  of  Figure  4-8  there 
are  5x4=20  pairs  (vj,v^)  initially  possible  and  of  these  only  10  agree  in  type 
of  tip  and  tail. 

The  mathematical  development  for  forming  a  quad  (A0,  Ax,  Ay,  As)  as  a 
function  of  v^  and  is  given  in  Figure  4-9.  Table  4-3  shows  computer  out¬ 

put  for  the  example  shown  in  Figure  4-8.  There  are  10  quads  produced  and  a 
cluster  of  size  3  is  evident  in  the  neighborhood  of  the  best  registration 
transformation  (A0=5.1O,  Ax=-75,  Ay=88,  As=0.5).  In  real  world  cases  there 
would  be  hundreds  of  quads  overall  and  a  few  dozen  in  a  cluster. 
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w . •  tx ■ i  - j^A^nr****  o*^*-*- .  • 


Map  Vectors 


(L, 170,220) 
(T, 100, 100) 
(L, 200,100) 
(L, 260, 70) 
(T, 150, 125) 


(X, 100, 200) 
(Y, 40, 150) 
(X, 220, 170) 
(X,  40,  70) 
(Y, 150, 50) 


Image  Vectors 

a.545,400)  -  (X, 200, 120) 
(T, 260, 240)  -  (Y, 100,245) 
(T, 140, 380)  -  (Y, 300, 380) 
(L, 420, 370)  -  (X, 360, 500) 


Figure  4-8.  a)  Example  set  of  map  vectors,  and 

b)  set  of  vectors  representing  an  image  to  be  registered 
to  the  map  by  RS&T  transformation. 


-  '--T~Tnwr  ~ru.. 


Assuming  that  vector  corresponds  to  vector  Vj  transformation 
parameters  (A0,Ax,Ay ,As)  are  gotten  as 

A0  -  0  -  0i 

As  =  length  of  /length  of  xk 

Ax  =  As  sinA0  —  As  A^^  cosAO  +  Cx 

Ay  =  -As  Ax  sinA0  -  As  A^.  cosA©  +  Cy 

The  resulting  registration  transformation  in  homogeneous  coordinates 
is 

As  cosA©  As  sinAG  0 

[u,v,l]  =  [x,y,l]  "As  sinA0  A s  cosA©  0 

Ax  Ay  1 

where  (x,y)  is  an  image  point  and  (u,v)  is  the  corresponding  map  point. 


Figure  4-9.  Mathematical  derivation  of  RS&T  transformation  parameters 
from  a  pair  of  vectors  and  V.  assumed  to  be  correspond¬ 
ing.  J 
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4. 2. 2.1  Example  wich  Scale 

The  following  experiment  demonstrates  the  utility  of  the  proposed  technique 
for  regions  where  cultural  activity  creates  features  such  as  straight  edges 
or  networks  of  lineals.  The  results  reported  are  typical  of  many  similar  ex¬ 
periments  . 

Point  features  were  identified  by  eye  on  an  aerial  image  from  the  mid- 
western  U.S.A.  using  two  different  measuring  devices  and  two  different  orien¬ 
tations.  Figure  4-10  shows  sample  imagery  while  Table  4-4  contains  .the  coor¬ 
dinates  of  the  selected  feature  points.  Points  are  labeled  according  to  the 
type  of  road  intersection  -  'L',  'T' ,  'X',  or  'Y'  -  which  they  represent  or 
are  labeled  'A'  indicating  an  arbitrary  point  on  the  road. 

50  vectors  were  chosen  to  model  the  map  of  Figure  4-10  while  64  vectors 
represented  the  image.  Tables  4-5  and  4-6  show  some  of  the  resulting 
quads  formed  by  matching  vectors  from  the  map  with  vectors  from  the 
image.  Note  that  stage  coordinates  have  been  divided  by  10  for  format  con¬ 
venience.  Table  4-5  shows  10  vectors  in  stage  coordinates  which  have  the  same 
type  of  tip  and  tail  as  the  vector  (4365,2747)  -  (2498,3641)  in  the  photo  model. 
Only  the  first  of  these  matches  is  correct  and  hence  only  the  first  quad  (5.39, 
-1020,4040.  5.01)  contributes  to  the  ultimate  cluster.  Of  the  50  x  64  =  3200 
pairs  (\K  ,v^)  possible  only  790  nroduce  quads  in  cluster  space  after  the  check 
on  tip  and  tail  type.  Table  4-6  shows  that  32  of  the  790  quads  form  a  cluster 
near  the  parameter  set  (AG=5.38,  Ax=-1 000.,  Ay=4000. ,  As=5.0).  Using  Table  4-6 
the  reader  can  verify  that  30  of  these  32  quads  represent  correct  vector  matches. 
For  instance,  the  first  quad  represents  the  matching  of  vectors  from  point  #2 
to  point  it  1  in  the  two  different  coordinate  systems.  The  asterisks  mark  the 
two  incorrect  matches  which  in  fact  are  outliers  in  the  set  of  32  quads.  A 
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Figure  4-10.  Aerial  photo  over  midwestern  U.S.A.  referred 
to  as  "4621".  Coordinates  of  labeled  points 
are  given  in  Table  4-4. 


Table  4-4 


Digitization  of  points  on  image  4621  using  Talos  digitizer 
(0.001  inch  resolution)  and  scanning  stage  (0.0005  inch 
resolution) . 


IMAGE  PT  # 

TALOS  * 

(photo) 

X  Y 

STAGE  * 
(transparency) 

X  Y 

DESCRIPTION  ** 

1 

2498 

3641 

5000 

5000 

X 

2 

4365 

2747 

8722 

6803 

T 

3 

5583 

2191 

11120 

8052 

X 

4 

3044 

4744 

3913 

7203 

T 

5 

3342 

5357 

3314 

8405 

T 

6 

3611 

5908 

2814 

9511 

T 

7 

1980 

6678 

- 

- 

T 

8 

3975 

5750 

3508 

9946 

T 

9 

5516 

5073 

6563 

11448 

T 

10 

7130 

3514 

10964 

12108 

Y 

11 

7430 

4091 

10367 

13304 

L 

12 

7723 

4700 

9866 

14600 

L 

13 

8642 

4258 

- 

- 

X 

14 

6952 

3223 

11209 

11461 

A 

15 

6942 

3613 

10565 

12000 

L 

16 

7233 

4181 

10010 

13101 

T 

17 

4861 

7575 

1710 

13600 

X 

18 

5980 

7667 

- 

- 

Y 

19 

6444 

7480 

- 

- 

T 

20 

9471 

6190 

- 

- 

T 

21 

5013 

7886 

1399 

14295 

A 

22 

4480 

7729 

- 

- 

T 

23 

7346 

9632 

- 

- 

A 

24 

9282 

5763 

- 

- 

A 

25 

9700 

6104 

- 

- 

Y 

26 

5279 

7999 

1513 

14804 

A 

27 

4010 

1991 

9388 

5480 

A 

28 

5480 

2552 

- 

- 

L 

29 

4808 

2527 

— 

A 

*  The  photo  in  Figure 

4-10  was 

rotated  roughly  50  degrees 

clockwise  when 

mounted  on  the  digitizer.  The  corresponding  transparency  was  not 
rotated  when  mounted  on  the  stage. 

**  L,X,T,  and  Y  indicate  type  of  intersection,  while  A  indicates  an 
arbitrary  point  on  a  straight  road  segment. 
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Table  4-5. 


First  10  quads  in  set  of  790  quads  produced  by 
comparing  50  map  vectors  with  64  image  vectors. 
(Stage  coordinates  of  Table  4-4  divided  by  10 
for  format  convenience). 


MAP  VECTOR  IMAGE  VECTOR  TRANSFORMATION 

VEC  TAIL  HEAD  VEC  TAIL  HEAD 


•  TYP(X,Y)  TYP(X,Y)  «  TYP(X,Y) 

t  3(4345,2747)  4(2498,3441)  51  3(  872,  480) 
1  3(4365,2747)  4(2498,3641)  52  3<  391,  7 20) 
t  3(4365,2747)  4(2498,3641)  S3  3(  331,  840) 
1  3(4365,2747)  4(2498,3641)  54  3(  281,  951) 
1  3(4365,2747)  4(2498,3641)  55  3(  350,  994) 
1  3(4365,2747)  4(2498,3441)  54  3<  654,1144) 
1  3(4365,2747)  4(2498,3641)  61  3(1001,1310) 
1  3(4365,2747)  4(2498,3641)  62  3(  872,  680) 
I  3(4345,2747)  4(2498,3641)  67  3(  872,  630) 
1  3(4365,2747)  4(2498,3441)  48  3(  391,  720) 


TYP(X, Y) 

THETA 

SCALE 

DELX 

DELY 

4  ( 

500, 

500) 

.539*01 

.501+01 

-.102*04 

.404+04 

4( 

500, 

500) 

.381*01 

.843+01 

.322*04 

.956*04 

4( 

500, 

500) 

.380+01 

.545*01 

.297*04 

.747*04 

4( 

500, 

500) 

.381+01 

.413*01 

.283*04 

.654*04 

4  ( 

500, 

50) 

.397*01 

.401*01 

.237*04 

.647+04 

4  ( 

500, 

500) 

.450*01 

.312*01 

.129*04 

.549*04 

4( 

500, 

500) 

.482*01 

.217*01 

.130*04 

.461+04 

4(1112, 

805) 

.221*01 

.745*01 

.125+05 

.534+03 

4( 

171, 

1360) 

.324*00 

.212*01 

.307*04 

.793+03 

4(1112, 

805) 

.258+01 

.285*01 

.640*04 

.389+04 

6« 


.»  *7  #***?**>& ujtm firf--- 


.-a’ 


32  quads  contained  in  the  bin  0£  [5. 0, 5.  7  ] ,  Se[4. 8,5.2] , 
Ax  e  [-1100,-900] ,  and  Ay  £ [3000 ,5000] . 


MAP  VECTOR 


IMAGE  VECTOR 


TRANSFORMATION 


VEC  TAIL 
#  TYP (X,Y) 


HEAD  VEC 

TYP(X,Y)  // 


TAIL 
TYP (X, Y) 


HEAD 

TYP(X, Y)  THETA  SCALE  DELX 


1 

2 

3 

4 
4 
6 
6 

7 

8 
9 

10 

11 

12 

19 

20 
21 
22 

24 

25 

29 

30 

31 
33 

35 

36 

37 

38 

39 

46 

47 

48 
50 


3(4365,2747) 

3(3044,4744) 

3(3342,5357) 

3(3611,5908) 

3(3611,5908) 

3(3975,5750) 

3(3975,5750) 

3(5516,5073) 

4(2498,3641) 

1(7430,4091) 

1(7723,4700) 

1(6942,3613) 

3(7233,4181) 

3(4365,2747) 

3(4365,2747) 

1(7430,4091) 

1(7723,4700) 

1(6942,3613) 

3(4365,2747) 

3(3044,4744) 

3(3342,5357) 

3(3611,5908) 

3(3975,5750) 

4(5583,2191) 

1(^430,4091) 

1(7723,4700) 

1(6942,3613) 

3(7233,4181) 

3(3044,4744) 

1(7430,4091) 

1(7723,4700) 

1(6942,3613) 


4(2498,3641) 

4(2498,3641) 

4(2498,3641) 

4(2498,3641 

4(2498,3641) 

4(2498,3641) 

4(2498,3641) 

4(2498,3641) 

5(7130,3514) 

4(2498,3641) 

4(2498,3641) 

4(2498,3641) 

4(2498,3641) 

4(5583,2191) 

5(7130,3514) 

3(4365,2747) 

3(4365,2747) 

3(4365,2747) 

4(4861,7575) 

4(5583,2191) 

4(5583,2191) 

4(5583,2191) 

4(5583,2191) 

5(7130,3514) 

4(5583,2191) 

4(5583,2191) 

4(5583,2191) 

4(5583,2191) 

5(7130,3514) 

3(3044,4744) 

3(3044,4744) 

3(3044,4744) 


51  3 (  872,  680) 

52  3(  391,  720) 

53  3(  331,  840) 

54  3 (  281,  951) 

55  3(  350,  994) 

54  3 (  281,  951) 

55  3 (  350,  994) 

56  3(  656,1144) 

57  4(  500,  500) 

58  1(1036,1330) 

59  1(  986,1460) 

60  1(1056,1200) 

61  3(1001,1310) 

62  3(  872,  680) 

63  3(  872,  680) 

64  1(1036,1330) 

65  1(  986,1460) 

66  1(1056,1200) 

67  3(  872,  680) 

68  3(  391,  720) 

69  3 (  331,  840) 

70  3(  281,  951) 

71  3(  350,  994) 

73  4(1112,  805) 

74  1(1036,1330) 

75  1(  986,1460) 

76  1(1056,1200) 

77  3(1001,1310) 

78  3(  391,  720) 

79  1(1036,1330) 

80  1(  986,1460) 

81  1(1056,1200) 


4(  500,  500) 
4(  500,  500) 
4(  500,  500) 
4(  500,  500) 
4 (  500,  500) 
4(  500,  500) 
4(  500,  500) 
4(  500,  500) 
5(1096,1210) 
4(  500,  500) 
4(  500,  500) 
4  (  500,  500) 
4(  500,  500) 
4(1112,  805} 
5(1096,1210) 
3(  872,  680) 
3 (  872,  680) 
3(  872,  680) 
4(  171,1360) 
4(1112,  805) 
4(1112,  805) 
4(1112,  805) 
4(1112,  805) 
5(1096,1210) 
4(1112,  805) 
4(1112,  805) 
4(1112,  805) 
4(1112,  805) 
5(1096,1210) 
3  (  391,  720) 
3  (  391,  720) 
3  (  391,  720) 


.539+01  .501+01-. 102+04 
.536+01  .501+01-. 101+04 
.536+01  .504+01-. 103+04 
.537+01  .504+01-. 104+04 
.553+01  .489+01-. 959+03 
.522+01  .514+01-. 994+03 
.538+01  .499+01-. 100+04 
.539+01  .504+01-. 105+04 
.538+01  .500+01-. 101+04 
.538+01  .501+01-. 102+04 
.538+01  .495+01-. 982+03 
.538+01  .497+01-. 992+03 
.538+01  .500+01-. 102+04 
.537+01  .495+01-. 941+03 
.538+01  .499+01-. 994+03 
.537+01  .499+01-. 987+03 
.538+01  .493+01-. 933+03 
.538+01  .493+01-. 923+03 
.538+01  .497+01-. 972+03 
.538+01  .496+01-. 962+03 
.537+01  .496+01-. 956+03 
.537+01  .499+01-. 992+03 
.538+01  .497+01-. 986+03 
.538+01  .502+01-. 105+04 
.537+01  .500+01-. 991+03 
.539+01  .494+01-. 959+03 
.538+01  .493+01-. 927+03 
.537+01  .500+01-. 101+04 
.538+01  .497+01-. 967+03 
.538+01  .499+01-. 991+03 
.538+01  .493+01-. 935+03 
.538+01  .495+01-. 955+03 


DELY 

.  404+04 
.412+04 
.411+04 
. 408+04 
.352+04 
. 464+04 
. 406+04 
.401+04 
. 404+04 
.407+04 
. 405+04 
.  406+04 
. 406+04 
. 408+04 
.405+04 
.410+04 
.402+04 
. 406+04 
. 406+04 
.407+04 
.410+04 
.410+04 
. 406+04 
. 407+04 
.414+04 
. 400+04 
. 404+04 
.410+04 
.404+04 
.  406+04 
. 406+04 
. 408+04 


4-D  clustering  routine  would  need  to  be  implemented  to  find  the  significant 
clusters. 


4.2.3  Extending  LNK  Registration  to  3  Dimensions 

So  far  the  LNK  registration  procedures  have  ignored  the  3D  aspects  of 
the  objects.  This  section  presents  a  method  to  extend  the  2D  registration 
to  3D  under  certain  constraints.  The  2D  procedure  found  RS&T  transforma¬ 
tions  mapping  map  structure  onto  image  structure.  Experiments  showed  that 
the  RS&T  assumption  was  feasible  in  cases  where  variation  in  the  3rd  dimen¬ 
sion  was  relatively  small.  Good  approximate  RS&T  registration  transforma¬ 
tions  were  obtained  automatically  for  photo/map  pairs  even  when  there  was 
some  relief  in  the  terrain.  There  are  many  cases  where  RS&T  transformations 
are  inadequate,  such  as  in  low  altitude  aerial  imaging  and  in  acquisition  of 
solid  objects  by  a  robot  vision  system.  In  these  cases  projective  transforma¬ 
tions  must  be  used. 

In  the  general  case  6  parameters  are  necessary  to  specify  imaging  in  a 
3D  world  [Duda  and  Hart  1973].  In  this  section,  we  consider  a  constrained 
imaging  system  with  only  3  free  parameters  as  shown  in  Figure  4-11.  A  front 
image  plane  is  used  with  reference  system  origin  (x=0,y=0,z=0)  at  the  image 
center.  The  camera  has  known  focal  length  f  and  looks  vertically  down  at  a 
scene  with  base  distance  yQ  from  the  image  plane.  There  are  only  3  unknown 
parameters  to  discover;  the  angle  0  at  which  the  object  lies  on  the  base 
plane  and  the  amount  of  translation  (xQ, z0)  of  the  object  origin  from  the 
point  where  the  camera  axis  pierces  the  plane. 

There  is  some  justification  for  this  assumption  in  the  aerial  imaging 
case.  First,  f  is  usually  known.  Second,  it  is  possible  to  get  a  good  ap¬ 
proximation  for  altitude  y0  and  to  achieve  a  nearly  vertical  camera  axis. 

These  approximations  would  perhaps  be  good  enough  to  correctly  detect  an 
approximate  3  parameter  view  which  could  be  used  as  an  initial  approximation 


Figure  4-11.  Three  dimensional  object  viewed  by  camera  with 

known  attitude.  Knowns  are  f  and  YQ  and  a  model 
of  the  object.  Unknowns  are  the  object  orienta¬ 
tion  parameters  Q,  X0,  and  ZQ. 
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to  a  full  6  parameter  view. 

We  therefore  proceed  as  follows.  A  3D  terrain  model  is  specified  in 
local  coordinates  as  in  Figure  4-12.  The  model  would  contain  the  location 
and  description  of  all  significant  features  such  as  edges,  corners,  inter¬ 
sections,  water  bodies,  etc. 

The  acquisition  problem  is  defined  as  discovering  (computing)  the  orien¬ 
tation  parameters  (9,x0,z0)  from  the  image  structure  and  the  known  camera 
parameters  f  and  yQ.  A  few  definitions  are  appropriate  before  proceeding. 

camera  parameters  -  parameters  that  fix  the  imaging  system  over  the 

base  plane,  i.e.  f  and  yQ,  and  define  the  global 
coordinate  system. 

orientation  parameters  -  parameters  that  fix  the  object  (or  object  model) 

in  the  global  coordinate  system  which  are  0,  x0, 
z0 . 

viewing  parameters  -  the  combined  camera  and  orientation  parameters  f,y0, 

0,xo,zo  which  allow  a  specific  image  to  be  created 
(from  the  terrain  model). 

Here  we  assume  that  we  know  the  attitude  of  the  camera  f,y0  and  the  ori¬ 
entation  0,xQ,zo  of  the  object.  We  develop  computational  formulas  for  image 
point  (x^,y^,z  )  corresponding  to  point  (xm»ym»zm)  the  map. 

Let  (x,y,z)  be  the  global  coordinates  of  point  (xm,ym,zm)  under  map 
orientation  (0,xo,zo).  Then  we  have 

(1)  x  =  x_  cos  0  -  z  sin  0  +  x„ 
m  m  o 

(3D  Map  to 


(2)  y  -  y0  -  ym 

(3)  z  =  xm  sin  0  +  zm  cos  0  +  zQ 


3D  Global) 


The  global  coordinates  (x,y,z)  then  produce  image  coordinates  (x^.y-pZ^) 
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via  the  direct  perspective  transformation. 


(4)  x  =  fx  =  f(*m  cos  0  ~  zm  sin  2  +  xo} 

‘ 1  f  +  y  f  +  y0  -  ym 

(5)  yr  =  o 

(6)  z  =  fz  -  f(xm  sin  9  +  zm  cos  6  +  zQ) 

1  f  +  y  f  +  y0  -  ym 

Here  we  assume  that  a  given  vector  (Ax,^  Aym,  Azm)  -  (Bx^,  Bym,  B^)  in 
the  map  corresponds  to  a  given  vector  (Ax’,  Az')  -  (Bx' ,  Bz')  in  the  image. 

We  develop  computational  formulas  for  determining  map  orientation  parameters 
(0,xo,yo)  from  this  correspondence. 

Rearranging  the  imaging  equations  (4)  and  (6)  from  above  we  have  the 
following. 


(7) 

XI  ^  f  *" 

+  y„> 

=  Xj,  cos 

6  - 

zm  sin  9  +  xo 

(8) 

N 

M 

1 

3 

+ 

•< 

o 

=  Xjj;  sin 

0  + 

zm  cos  9  +  z0 

f 

Since  (7)  must  hold  for  both  points  A  and  B,  we  get  two  equations  from 
which  xQ  can  be  eliminated,  leaving  only  S  unknown. 

(9)  AxI(f  -  Aym  +  y0)  -  Bx-^f  -  Bym  +  y0)  =  (Aj^  -  Bx^cos  9  + 

f  f 

(Bzm  -  Azm)sin  9 
Equation  (9)  is  of  the  form 

c  =  d  cos  9  +  e  sin  9 

where  we  make  the  substitutions  w  =  sin  0  and  /l-w^  =  Cos  9. 

Thus  a  standard  quadratic  equation  in  w 

(10)  (e2  +  d2)w2  -  (2ce)w  +  (c2  -  d2)  =  0 

is  obtained  where  the  coefficients  are  gotten  as  follows: 


(3D  Global 

to 

2D  Image) 
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c  =  (Axx  (f  -  Aym  +  y0)  -  BxT ( f  -  Bym  yQ))/ 

d  =  Axm  “  Bxm 
e  =  Bzm  -  Azm  • 

Solving  the  quadratic  yields 

(11)  w  =  ce  +  d  Ve2  -  c2  +  d2 


(12)  0  »  sin-d  w. 

Knowing  6,  (7)  and  (8)  can  be  used  to  solve  for  xQ  and  yQ  using  either  A  or 
B  point  coordinates.  For  example, 

xq  -  Axr(f  -  Aym  +  y0)/f~Axm  cos  6  +  Azm  sin  0. 

Since  there  is  mathematical  ambiguity  in  6  from  (11)  two  parameter  sets 

(wi*x01,y0l)  and  ^w2,X02’y02^  Can  resu^c*  is  easY  to  check  for  correct 

parameters  using  (4)  and  (6)  and  the  2  known  pairs  of  corresponding  points 

(AI-V  and  (BI’V- 

There  are  two  significant  cases  to  note  where  the  computation  breaks 
down.  First  of  all,  the  discriminant  of  the  quadratic  can  be  negative  and 
hence  no  0  can  be  obtained.  This  will  happen  whenever  the  map  edge  cannot 
possibly  be  imaged  onto  the  image  edge.  Few  pairings  are  actually  possible 
due  to  the  fixed  scale  imposed  by  f  and  yQ.  Secondly,  whenever  the  map  edge 

is  vertical  both  d  and  e  above  are  zero  and  equations  (11)  and  (12)  cannot 

produce  0.  Physically  we  can  rotate  the  object  in  the  map  freely  about  that 
vertical  edge  without  altering  its  image  and  thus  we  should  not  expect  to  get 
-*  mathematically  either. 

Figure  4-13  shows  a  planar  section  containing  the  camera  axis,  the  ver- 
i map  vector,  and  hence  the  image  vector  as  well.  Clearly,  free  rotation 
••  -.-ct  about  the  axis  AB  will  not  change  this  picture.  It  is  also 

r ,  that  locational  parameters  x  and  z  are  completely  specified 
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by  A'B'  and  AB  provided  that  they  correctly  match  via  the  viewing  transforma¬ 
tion.  From  the  imaging  equations  (4)  and  (6)  applied  to  both  points  A  and  B 
we  get 

AXj. (f+yQ)  /  f  =  x  =  BXjCf+y)  /  f  or 

(23)  AxI(f+yQ)  =  BXjCf+y) 

and  similarly  for  the  z  coordinate  relations 

(24)  AzI(f+yo)  =  BzI(f+y). 

Conditions  (23)  and  (24)  must  hold  if  A'B'  is  to  possibly  match  with  AB. 

We  are  already  in  trouble  here  if  only  real  edge  segments  are  a-  .liable 
because  (23)  and  (24)  are  scaling  equations.  On  the  other  hand,  if  abstract 
vertical  map  edges  which  have  accurate  tips  and  tails  are  being  used,  it  would 
be  easy  to  consider  only  nonvertical  vectors  as  before.  All  that  need  be  done 
is  to  construct  nonvertical  vectors  by  mixing  the  tip  and  tail  points  of 
several  vectors.  For  this  reason,  special  treatment  of  vertical  vectors  can 
be  ignored  (justifiably)  in  the  computer  programs. 

As  an  example,  suppose  the  map  contained  the  object  shown  in  Fipure  4-12. 
Suppose  the  image  was  formed  by  viewing  the  object  with  aspect  parameters 
(f=l.yo“10)  and  orientation  parameters  (0=30° ,xQ=-2,zo=3) .  The  resulting 
image  is  shown  in  Figure  4-14,  along  with  the  coordinates  of  the  points.  Let 
the  object  in  the  map  and  image  be  represented  by  the  vectors  listed  in  Table 
4-7.  Matching  the  10x10  pairs  of  vectors  (v^v^)  yields  only  12  feasible  para¬ 
meter  sets,  10  of  which  form  a  cluster  about  the  correct  parameters  a  *  (0=0.525, 

xo=-2, zQ= 3 ) . 
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4.2.4  Clustering 

All  of  the  versions  of  LNK's  registration  technique  require  clustering  to 
be  done  in  a-space.  The  cluster  space  is  filled  with  points,  each  one  of  which 
represents  the  matching  of  one  image  element  to  one  map  element  on  the  basis  of 
local  features  alone.  A  cluster  of  points  in  this  space  represents  a  globally 
good  transformation  that  matches  many  image  elements  to  corresponding  map  ele¬ 
ments.  A  threshold  may  be  set  to  determine  the  strength  of  an  acceptable 
cluster;  for  instance,  we  might  demand  that  half  of  the  map  elements  be  matched 
by  image  elements.  A  unique  cluster  will  typically  result  but  there  may  be  no 
cluster  formed  in  the  case  where  there  is  poor  feature  detection  or  where  the 
image  and  map  really  do  not  match.  On  the  other  hand,  several  strong  clusters 
can  result  as  in  the  case  where  symmetry  produces  several  good  matching  possi¬ 
bilities.  LNK  has  used  two  different  clustering  techniques  -  hierarchical 
clustering  and  variable  resolution  clustering. 

In  the  hierarchical  technique,  clustering  for  a=(A0,Ax,Ay)  was  done,  first 
on  Q  alone,  and  then  in  (Ax, Ay)-space  given  a  fixed  A0.  This  technique  was 
useful  for  human  interaction  since  0-space  could  be  viewed  as  a  histogram  and 
(Ax,Ay)-space  as  a  scatter  plot. 

The  0-space  was  divided  into  360  one-degree  bins  and  the  rotation  A0^m  that 
rotated  the  image  edge  to  the  map  edge  was  entered.  After  smoothing,  the  top 
three  peaks  of  the  histogram  were  chosen  for  the  second  step.  For  each  of  the 
transformations  that  had  a  A0  in  one  of  the  three  peaks,  the  (Ax.Av)  in  that 
transformation  was  recorded  in  (Ax,Ay)-space.  The  resulting  scatter  plot  could 
then  be  searched  for  clusters. 

Intuitively,  clustering  can  be  done  by  moving  a  box  around  the  cluster 
space  to  see  if  there  is  some  position  at  which  an  acceptable  number  N  of  points 
are  inside  the  box.  If  so,  the  coordinates  (G,s,Ax,Ay)  of  that  position  mark  a 
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cluster  center  and  represent  a  good  registration  transformation.  A  computer 
implementation  can  be  specified  when  we  choose  the  size  of  the  box  and  the  set 
of  different  positions  (bins)  at  which  we  will  try  to  fill  it. 

The  ultimate  spread  of  a  cluster  is  related  to  the  precision  or  resolu¬ 
tion  of  the  imaging  system  and  the  feature  detectors.  Suppose  the  Hough  de¬ 
tector  is  used  to  detect  straight  edge  elements  and  the  precision  is  ±1  degree 
in  direction.  An  error  of  1  degree,  or  fr/180  radians,  in  0  could  result  in  an 
error  of  rTr/180  in  the  Ax  or  Ay  computed  from  0  as  in  Figure  4-9,  where  r  is 
the  distance  from  the  origin  to  the  liAe  formed  by  the  straight  edge.  If  r  is 
900  pixels  then  the  error  induced  in  Ax  and  Ay  could  be  15  pixels.  Thus  in  this 
case  the  box  (bin)  size  should  never  be  smaller  than  2  degrees  x  30  pixels  x  30 
pixels.  In  the  variable  resolution  technique  binning  implementation  grids  of 
10  x  10  x  10  boxes  were  used.  Each  point  (0,Ax,Ay)  was  distributed  to  the  ap¬ 
propriate  bin.  Assuming  that  the  ranges  of  0,Ax,  and  Ay  were  360  degrees,  2000 
pixels,  and  2000  pixels  respectively,  the  bin  size  of  level  one  would  be  36 
degrees  x  200  pixels  x  200  pixels.  If  any  of  the  1000  bins  at  level  one  received 
N  or  more  triples  binning  was  repeated  on  only  those  points  at  level  two,  and 
with  a  bin  size  of  3.6  degrees  x  20  pixels  x  20  pixels.  A  third  level  of  cluster¬ 
ing  was  seldom  warranted  and  in  the  case  at  hand  would  imply  a  box  size  of  0.36 
degrees  x  2  pixels  x  2  pixels. 

There  is  one  final  note  on  the  clustering  procedure.  Since  the  cluster 
space  is  quantized  by  the  bins,  overlapping  sets  of  bins  are  actually  used  so 
that  clusters  are  not  lost  by  being  split  by  bin  boundaries.  Thus  each  triple 
(Q,Ax,Ay)  is  actually  distributed  to  several  bins  rather  than  one. 

To  use  detected  straight  edges  instead  of  abstract  edges,  some  changes  need 
to  be  made  in  the  clustering  algorithm.  Real  edge  detectors  create  edges  with 
two  general  defects.  First  of  all,  only  small  segments  of  the  true  edge 
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are  usually  detected  primarily  because  only  a  fixed  size  detector  is  used. 


Second,  the  detectors  may  overshoot  the  true  termination  of  an  edge  at  a 
corner. 

In  order  to  deal  with  the  above  problems  when  real  edges  are  used,  the 
matching  of  image  vectors  to  map  vectors  must  be  made  to  be  sloppy.  Figure 
4-15  gives  the  details.  The  image  edge  element  must  be  allowed  to  slide  along 
the  potentially  matching  model  edge.  Mathematically  this  generates  an  infinite 
set  of  triples  (6, Ax, Ay)  constrained  to  be  on  a  line  segment  in  cluster  space. 
In  practice  the  endpoints  (G.Ax-j-.Ayj.)  and  (G.Ax^.Ay^)  of  the  line  segment,  are 
saved  for  clustering.  Whenever  points  (Q,Ax.,Ay^)  and  (Q,Ax9,Ay2)  are  in  sepa¬ 
rate  bins,  Doints  are  also  contributed  to  all  bins  in  between  them.  This  tech¬ 
nique  was  used  for  all  real  edge  vectors  in  the  experiments  reported  in  Section 
4.2. 
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Image  edge  AB  or  (Q^,^)  is  rotated  0,,,-^  onto  A'B'  or  (0m»ri) 
to  be  parallel  to  map  edge  CD  or  (0m, rm) .  Projection  of  A'C  onto 
A'R  has  constant  length  yielding  equation  relating  Ax  and  Ay 
(translation  parameters) 

Ax  cos  0  +  Ay  sin  0„  +  (r.-r  )  =  0. 
m  '  m  j.  m 

(Ax, Ay)  are  further  constrained  because  point  A'  should  translate 
no  further  than  point  C  and  point  Bf  should  translate  no  further 
than  point  D.  Some  tolerance  can  be  given  on  either  end  since 
edge  detectors  can  overshoot  corners. 


Figure  4-15.  Constraints  on  points  a  =  (G.Ax.Ay)  in  cluster 
space  derived  by  pairing  real  image  edge  segment 
A3  with  real  map  edge  segment  CD. 


4. 3  Region  Matching 

Symbolic  matching  of  images  consists  of  reducing  two  images  to  a  a  pair 
of  abstract  structures  and  comparing  these  structures.  The  methods  we  con¬ 
sider  here  take  the  regions  in  a  segmented  image  as  components  of  the  struc¬ 
ture  describing  the  image.  Relationships  between  regions  may  be  stored  ex¬ 
plicitly  in  a  regions  description,  or  they  may  represented  as  edges  in  a  graph. 
Region  descriptions  may  include  such  features  as  size,  shape  and  average  gray 
level.  All  matching  procedures  based  on  regions  rely  on  the  segmentation 
procedures  used  to  define  the  regions.  Poor  region  segmentation  can  cause 
serious  mismatching  of  regions.  The  choice  of  features  describing  regions  may 
also  be  critical  in  the  design  of  those  procedures. 

Most  region  matching  algorithms  do  not  provide  the  registration  trans¬ 
formation.  This  means  additional  processing  must  be  done,  after  the  region 
matching  procedures  are  through,  to  find  the  desired  transformation.  Since 
region  segmentation,  region  features  and  region  matches  can  all  be  computa¬ 
tionally  expensive  to  obtain,  region  matching  should  be  approached  with  cau¬ 
tion.  The  expense  might  be  justified  if  image  interpretation  is  desired 
where  region  descriptions  could  help  tremendously. 

The  output  of  a  region  matching  procedure  applied  to  two  images  is  a  cor¬ 
respondence  between  some  subset  of  the  set  of  regions  in  the  first  image  and  a 
subset  of  the  set  of  regions  in  the  second  image.  If  we  label  the  regions  in 
image  1  with  the  symbols  a^,  ...,  aN  and  the  regions  in  image  2  with  the  symbols 
b^,  ...,  b^,  then  a  matching  assigns  an  element  of  the  set  {b^,  ...,  bM,  c}  to 
each  a^,  i=l,  . . . ,  N.  If  bj  is  assigned  to  a^  the  corresponds  for  ai.  If  c 
is  assigned  to  a^  then  a^  doesn't  correspond  to  a  region  in  image  2.  We  have 
been  unable  to  find  algorithms  in  the  literature  for  determining  a  registration 
from  a  matching  of  regions. 
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Prior  to  constructing  such  algorithms,  we  must  have  a  criterion  for  the  quality 
of  a  registration  obtained  from  region  matching.  One  natural  criterion  is  the 
maximization  of  the  total  area  matched  under  a  registration.  If  region 
matches  with  b j ,  then  under  a  registration,  we  can  view  both  a^  and  bj  as  regions 
in  the  same  plane  so  their  matched  area  is  a^  H  bj . 

If  a  measure  of  matched  quality  based  on  feature  values  of  a ^  and  bj  is 
available,  then  the  matched  area  may  be  weighted  by  the  match  quality.  This 
criterion  could  be  used  with  a  hill  climbing  optimization  technique  and  an 
initial  registration  to  provide  a  search  technique  for  a  registration.  This 
procedure  could  be  very  costly  if  the  initial  registration  is  poor. 
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4.3.1  Region  Image  Matching  Using  Similarity  of  Region  Features 

Region  image  matching  can  be  done  by  finding  a  set  of  features  which  de¬ 
scribe  the  regions  and  then  pairing  regions  which  have  the  best  matching  set 
of  features.  It  is  desirable  that  the  decisions  take  into  account  the  adja¬ 
cency  information  of  the  regions,  so  that  adjacent  regions  in  the  image  match 
adjacent  regions  in  the  map. 

One  such  method  of  region  matching  yields  a  measure  of  similarity  between 
pairs  of  regions,  one  from  each  of  the  two  images  to  be  matched  [Price  and 
Reddy  1979].  This  method  makes  no  assumptions  about  the  relative  displace¬ 
ment  and  orientation  of  the  pictures.  The  steps  of  their  algorithm  are  as 
follows: 

1)  Segment  the  image. 

2)  To  each  region  i  assign  a  set  Vil*  Vin  where  denotes  the 

value  of  the  j th  feature  for  region  i.  These  numbers  may  describe 
features  such  as  shape,  size,  position,  spectral  values,  etc.  To 
each  pair  of  regions,  region  i  from  image  1  and  region  j  from  image 
2,  define  the  region  to  region  match  rating,  R^ ,  by: 

8u  -  -  £  |  vik  -  v  |  \  sk 

where  is  a  normalization  factor  for  the  k^  feature  and  is  a 
measure  of  importance  of  the  k1"^  feature.  Larger  values  of  R^j 
indicate  good  matches. 

3)  In  this  step  we  attempt  to  improve  the  accuracy  of  the  rating 

by  taking  into  account  adjacency  information.  To  each  region  R^  in 
image  1,  assign  the  region  Rj  in  image  2  which  maximizes  j .  Let 
denote  the  set  of  regions  in  image  1  which  are  neighbors  of  region 
R^  and  have  a  match  in  image  2.  We  say  two  regions  are  neighbors  if 
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they  have  a  common  boundary  point.  Let  Nj  denote  those  regions  in 


image  2  which  match  the  regions  in  N^  Let  Nj '  denote  those  regions 
in  N.  which  are  neighbors  of  Rj .  The  neighbor  feature  value  (NFV)  of 
is  defined  to  be  the  number  of  elements  in  N^  and  the  NFV  for  Rj 
is  the  number  of  elements  in  N^ ' .  Recompute  the  Rjj's  using  this 
additional  feature  and  assign  to  each  region  in  image  1  the  region 
in  image  2  which  matches  it  best. 

This  matching  algorithm  is  designed  to  be  invariant  under  rotation  and 
translation  of  the  images.  By  omission  of  size  features,  the  algorithm  can 
be  made  invariant  to  scale  change.  Dissimilar  image  matching  and  image  to  map 


matching  can  be  handled  by  this  method. 


4.3.2  Region  Adjacency  Graph 


One  method  of  symbolic  matching  is  to  represent  region  segmented  images 
I  or  maps  as  graphs  and  to  use  graph  theoretic  techniques  to  do  the  matching 

[Pavlidis  1977].  A  graph  is  a  set  whose  elements  are  called  nodes  and 

|  a  collection  of  pairs  from  {v^}  called  edges.  For  scene  matching,  each  node 

^  corresponds  to  a  region  and  an  edge  in  the  graph  corresponds  to  two  regions 

which  are  adjacent  in  the  image.  The  graph  thus  formed  is  called  the  region 

I  adjacency  graph  (RAG)  of  the  image.  Each  node  contains  a  set  of  labels  de¬ 

scribing  the  region  represented.  These  labels  give  information  about  shape, 

|  size,  average  brightness,  etc.  The  RAG  of  a  map  database  may  contain  addi- 

g  tional  information  in  the  node  labels  such  as  the  type  of  region,  e.g.  forest, 

fields . 

I  Denote  the  nodes  of  the  RAG  for  image  1  by  {nl}  and  the  nodes  of  the  RAG 

for  image  2  by  {m^}.  We  ultimately  wish  to  match  elements  of  with  elements 

|  of  .  Thus  we  must  define  a  notion  of  matching  between  the  labels  attached 

|  to  nodes  and  for  any  i  and  j.  The  notion  of  label  matching  is  problem 

dependent,  but  typical  matching  conditions  for  two  regions  would  be  to  have 
I  approximately  the  same  area,  shape  and  average  brightness.  Poor  segmentation 

may  necessitate  a  more  involved  definition  of  matching.  We  would  also  like  to 
|  require  the  condition:  if  node  matches  and  matches  Mp  and  and 

■  are  adjacent,  then  M.  and  M  should  be  adjacent. 

1  3  P 

The  problem  of  finding  the  largest  number  of  matching  nodes  is  called  the 
I  subgraph  isomorphism  problem  and  is,  in  general,  computationally  expensive. 

However  good  segmentation  and  well-chosen  labels  could  reduce  the  cost  con- 
|  siderably.  In  addition,  the  cost  could  be  reduced  by  ignoring  adjacencv  in- 

■  formation  or  by  finding  a,  not  necessarilj^  good,  approximate  registration  first 
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In  the  latter  case,  a  prerequisite  for  label  matching  is  that  regions  should 
be  in  approximately  the  same  location  under  the  registration. 

In  either  case,  an  approximate  solution  can  be  obtained  using  a  maximal 
network  flow  algorithm  (Frank  and  Frishch,  1971].  Under  either  of  the  above 
assumptions,  the  RAG's  for  the  two  images  have  no  edges.  To  apply  the  maximal 
flow  algorithm,  we  form  a  new  graph  G  with  vertices  {M^f  U  {sourcef  U  {sink} 

where  source  and  sink  are  two  new  nodes.  There  is  an  edge  from  source  to  each 
M^,  and  an  edge  from  each  to  the  sink.  For  each  i  and  j  there  is  an  edge 
from  to  if  node  can  be  matched  with  .  In  addition,  we  assign  the 
integer  one,  called  the  capacity  or  weight,  to  each  edge  of  the  graph.  Finally 
an  optimal  flow  algorithm  is  applied  to  the  graph.  The  output  of  this  program 
is  an  assignment  of  either  zero  or  one  to  each  edge  of  the  graph.  If  we  erase 
the  source  and  sink  and  all  edges  which  have  been  assigned  a  value  zero,  then 
the  remaining  graph  will  have  the  property  that  no  node  has  more  than  one  edge 
incident  to  it.  If  we  erase  nodes  having  no  edges,  then  the  remaining  edges 
define  a  one-to-one  correspondence  between  a  subset  of  the  regions  of  one  image 
with  a  subset  of  the  regions  of  the  other  image.  This  correspondence  has  the 
property  that  it  maximizes  the  number  of  regions  matched. 

The  above  method  assumed  that  each  region  in  an  image  could  correspond 
to  at  most  one  region  in  the  map  or  other  image.  Under  real-world  conditions, 
this  assumption  is  unlikely  to  hold  as  it  would  be  quite  possible,  as  an  example, 
for  a  forest  which  is  one  region  in  a  map,  to  be  segmented  into  several  regions 
in  the  image.  The  above  algorithm  could  handle  this  possibility  by  changing 
the  initial  weights  of  the  edges,  connecting  the  sink  with  nodes  ,  to  a 
number  greater  than  one.  If  the  {m^}  are  the  map  regions  and  |n^}  are  the 
image  regions,  this  will  allow  a  many-to-one  relation  from  the  map  to  the  Image. 
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The  RAG  is  designed  to  be  invariant  under  rotation  and  translation  of 
the  images.  If  features  used  are  based  on  shape  and  not  size,  then  the  RAG 
is  also  invariant  under  scale  change.  The  RAG  can  be  used  for  matching  an 
image  to  a  map  or  for  matching  images  from  dissimilar  sensors. 


4.4  Comparison  of  Registration  Procedures 


Not  all  of  the  methods  could  easily  handle  the  full  RS&T  transforma¬ 
tions.  While  the  LNK  registration  procedure  can  easily  handle  the  full  RST 
transformation,  most  of  the  other  algorithms  could  not.  Of  the  correlation- 
based  methods,  only  the  invariant  moment  correlation  is  invariant  under  RST 
transformations.  It  is  possible  to  extend  the  other  correlation-based  meth¬ 
ods  to  include  rotation  by  repeating  the  correlation  while  rotating  the  image 
with  respect  to  the  map.  Except  for  the  hierarchical  technique,  where  the 
rotation  could  be  done  at  low  resolutions,  this  extension  would  be  extremely 
time  consuming.  The  region  matching  algorithms,  with  careful  selection  of 
region  features,  are  also  invariant  under  RST  transformations. 

The  space  requirements  for  storing  the  map  varies  greatly,  even  within 
one  registration  technique,  depending  on  the  features  being  used.  With  cor¬ 
relation,  if  a  whole  grayscale  Image  is  used  as  the  map,  then  there  is  no 
data  compression.  Whereas  if  edges  are  used  for  correlation-based  techniques 
or  LNK's  registration  procedure,  then  only  the  chain  codes  (if  the  edges  are 
not  straight)  or  the  endpoints  (if  the  edges  are  straight)  need  be  stored. 
This  could  lead  to  a  fair  amount  of  data  compression.  If  only  point  features 
are  used  with  correlation  or  LNK's  procedure,  then  a  fair  amount  of  data  com¬ 
pression  is  again  possible. 

The  space  requirements  for  a  map  using  a  region  matching  procedure  are 
highly  variable.  It  largely  depends  on  the  number  and  type  of  features  used 
to  characterize  the  regions.  What  is  needed  to  be  stored  is  the  adjacency 
relationships  and  the  features  characterizing  each  region.  As  long  as  none 
of  the  features  are  space  consuming,  such  as  the  chain  code  of  the  boundary, 
then  the  data  compression  could  be  large. 
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In  all  of  the  techniques,  the  amount  of  data  compression  must  be  bal¬ 
anced  against  the  need  for  accuracy.  It  is  better  to  allow  the  map  to  have  a 
large  number  of  features  and  let  the  feature  extraction  on  the  image  be  less 
extensive  but  reliable.  This  means  that  any  image  feature,  which  should  be 
reliable,  would  have  a  corresponding  map  feature  but  not  necessarily  vice 
versa. 

The  relative  speeds  of  the  registration  techniques  are  highly  variable. 

It  is  expected  that  LNK's  registration  procedure  would  perform  the  fastest; 
straight  correlation  would  be  very  slow.  Sequential  and  hierarchical  correl¬ 
ation,  while  three  orders  of  magnitude  faster  [Hall  1979]  than  straight  cor¬ 
relation,  would  be  slower  than  LNK's  procedure.  Since  region  matching 
methods  require  image  segmentation,  they  are  potentially  very  slow.  In  addi¬ 
tion,  since  more  processing  must  be  done  to  find  the  registration  transforma¬ 
tion,  region  matching  becomes  a  very  slow  process. 

The  robustness  of  the  registration  techniques  is  an  important  feature. 

In  all  cases,  the  robustness  of  any  registration  technique  is  going  to  depend 
on  the  robustness  of  its  corresponding  feature  detector.  It  is  difficult  to 
say  much  about  the  robustness  of  the  procedures  without  experimentation,  but  a 
few  things  are  known.  It  is  known  that  correlation  techniques  using  gray  scale 
images  are  highly  susceptible  to  noise.  Correlation  techniques  using  edges 
which  are  not  fattened  also  have  difficulties. 

LNK  has  performed  enough  experimentation  with  the  LNK  registration  pro¬ 
cedure  to  know  that  with  mediocre  feature  detection,  the  registration  technique 
performed  well.  The  region  matching  procedures  are  very  heavily  dependent  upon 
the  region  segmentation.  Poor  region  segmentation  would  make  region  matching 
almost  impossible  for  most  cases. 

It  might  be  justifiable  to  use  a  more  time  consuming  or  more  space  con- 
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suming  registration  procedure  if  the  extra  time  or  space  produces  extra 
information  which  can  later  be  used  for  scene  analysis.  For  example,  while 
finding  and  classifying  intersections  is  more  time  consuming  than  finding 
just  edges,  it  can  help  in  later  identification  of  roads  and  buildings. 
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4.5  Image  Disparity  Determination 

Once  a  global  registration  has  been  determined,  it  usually  is  necessary 
to  modify  the  transformation  to  account  for  local  distortion.  The  distortion 
could  be  caused  by  motion,  temporal  scene  changes,  imperfect  image  rectifica¬ 
tion,  or  local  altitude  variations.  The  disparity  between  two  registered 
images  is  the  small  local  differences  in  corresponding  pixel  locations  caused 
by  these  distortions.  The  analysis  of  image  disparity  in  similar  scenes  is  a 
key  element  in  the  computer  analysis  of  motion  as  well  as  in  the  construction 
of  depth  information  from  stereograms. 

Barnard  and  Thompson  [1980]  give  an  algorithm  which  takes  as  input,  two 
images  and  an  approximate  global  registration,  and  outputs  a  set  of  point  fea¬ 
tures  from  each  image  and  a  1-1  map  identifying  the  correspondence  between 
these  two  sets.  This  correspondence  can  be  used  to  define  a  piecewise  linear 
transformation  between  the  images  which  refines  the  global  registration.  Three 
properties  of  images  play  a  critical  role  in  determining  the  success  of  their 
approach. 

(1)  The  measure  of  distinctness  of  points. 

(2)  The  measure  of  similarity  of  such  points  taken  from  the  two  images. 

(3)  The  measure  of  how  well  a  point  match  accords  with  matches  of  neigh¬ 
boring  points. 

Point  features  lie  in  highly  variable  areas  in  which  the  variation  is  large 
in  all  directions  from  the  point.  An  interest  operator  given  by  Moravec  [1977] 
defines  a  point  to  be  of  high  interest,  if  there  is  a  large  variance  of  the 
gray  levels  along  the  horizontal,  vertical  and  diagonal  directions  within  the 
5x5  neighborhood  about  the  point.  The  interest  operator  tentatively  assigns 
the  minimum  of  these  variances  as  the  interest  value  of  the  point.  Finally 
local  maxima  with  interest  values  above  some  threshold  are  selected  as  the 
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true  points  of  interest. 

The  approximate  global  registration  is  applied  to  the  two  images.  A  re¬ 
laxation  procedure  is  then  applied  to  match  points  in  the  two  images.  Each 
interest  point,  a,  in  image  1  is  assigned  a  set  of  labels  and  a  corresponding 
set  of  disparity  vectors.  The  labels  are  the  points  (b^}  in  image  2  lying  in 
a  small  neighborhood  about  a’s  transformed  location.  The  disparity  vector 
corresponding  to  label  b^  is  the  vector  from  point  a  to  point  b^.  The  merit 
of  a  label  b^  is  based  on  the  similarity  of  the  5x5  windows  surrounding  a  and 
b^.  An  updating  rule  is  then  applied  to  the  merits  as  follows.  The  merit  of 
a  label  is  increased  when  nearby  points  in  image  1  have  labels  of  high  merit 
with  disparity  vector  similar  to  that  of  b^.  Barnard  and  Thompson  considered 
ten  iterations  of  this  procedure  sufficient. 

This  procedure  is  not  restricted  to  points  derived  using  Moravec's  interest 
operator.  Any  of  the  point  features  described  in  section  3.3  could  be  used.  In 
particular,  the  Moravec  interest  operator  is  unlikely  to  work  well  when  register¬ 
ing  images  from  different  sensors,  since  it  is  dependent  upon  gray  scale.  How¬ 
ever,  problems  involving  registering  images  to  maps  may  be  solved  using  this 
method. 

In  the  case  of  matching  images  from  different  sensors,  the  notion  of  similar 
neighborhoods  for  corresponding  points  would  have  to  be  redefined  since  gray 
levels  may  not  match  well.  Similarity  of  neighborhood  texture  might  provide  a 
feasible  alternative  to  grey  level  matching.  In  all  cases,  if  the  feature  points 
are  too  sparse  in  the  images,  the  procedure  for  updating  the  merits  will  require 
considerable  revision  since  small  neighborhoods  about  the  feature  points  may  not 
contain  other  feature  points. 
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5.  Conclusions  and  Recommendations 

The  preceding  sections  report  on  L.N.K.'s  study  of  the  problem  of 
registration  of  images  with  maps  or  dissimilar  images.  The  registration  pro¬ 
cess  is  divided  into  four  steps: 

1.  Determining  and  extracting  the  appropriate  features. 

2.  Computing  the  parameters  of  an  approximate  global  registration 
transformation. 

3.  Finding  the  disparity  for  a  subset  of  points  in  the  images. 

4.  Determining  a  global  nonlinear  transformation  for  the  entire  image. 

Three  classes  of  features  and  their  detectors  are  presented  in  this  re¬ 
port:  edge  features,  point  features,  and  region  features.  Various  methods  for 
registration  are  described  in  the  preceding  sections  including  several  correla¬ 
tion  procedures,  the  L.N.K.  registration  technique  and  its  extension,  and  some 
region  matching  procedures.  The  registration  techniques  are  briefly  compared 
and  contrasted.  Their  ability  to  handle  the  full  RS&T  transformations  and  their 
time  and  space  requirements  are  examined. 


5. 1  Conclusions 


Based  on  L.N.K.'s  experience  with  feature  detection  and  registration  and 
the  reported  success  of  various  algorithms  presented  in  the  literature,  the 

conclusions  of  this  study  are  the  following: 

.  Edge  features  or  edge-based  point  features,  such  as  intersections  and 
high  curvature  points,  are  the  most  promising  features  for  registra¬ 
tion. 

.  Abstract  edges  and  triangles  formed  from  point  features  are  very  use¬ 
ful  for  registration.  They  can  provide  a  more  accurate  initial  trans¬ 
formation  for  little  additional  feature  detection  effort.  In  particular 
the  clustering  in  L.N.K.'s  registration  procedure  will  be  performed 
faster  using  abstract  edges  than  real  edges. 

.  Only  two  of  the  registration  techniques  studied  are  worthy  candidates 
for  further  investigation.  L.N.K.'s  registration  procedure  can  pro¬ 
vide  a  full  RS&T  transformation  and  still  work  well  even  with  only  fair 
feature  detection.  The  other  promising  technique  is  a  combined  hierar¬ 
chical-sequential  correlation  method  using  edge  images. 

.  Map-guided  registration  can  aid  in  the  registration  process  by  provid¬ 
ing  as  part  of  its  content  what  feature  to  use,  the  particular  feature 
detector  and  appropriate  window  size  and  threshold  information,  and 
which  registration  method  to  use.  (For  example,  if  the  map  contains  a 
scene  with  man-made  structures  then  intersections  would  be  a  good  fea¬ 
ture,  whereas  in  scenes  with  few  man-made  structures  high  curvature 
points  could  be  used.)  The  use  of  maps  in  guiding  the  decision  making 
enables  the  process  to  be  tailored  to  the  given  situation. 

.  The  determination  of  the  disparity  for  a  subset  of  points  in  the  image 
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is  an  important  open  problem.  (One  method  is  described  in  this  report  but 
much  more  work  on  this  subject  needs  to  be  done.) 


5.2  Recommendations 


Based  on  this  study  the  following  theoretical  and  experimental  investiga¬ 
tions  are  proposed. 

I.  L.N.K.  recommends  that  the  following  feature  detection  algorithms  be  im¬ 
plemented  on  the  DIAL  system  and  tested: 

a.  Hough  edge  detector 

b.  Laplacian  (of  the  Gaussian  function)  edge  detector 

c.  Linking  of  edges  algorithm 

d.  Intersection  detection  algorithm 

e.  High  curvature  points  algorithm 

f.  Moravec  interest  operator  algorithm 

II.  The  accuracy  and  robustness  of  a  registration  procedure  are  dependent 
upon  the  corresponding  feature  detectors.  Various  facets  of  this  depen¬ 
dency  need  to  be  investigated.  Two  of  the  factors  are  the  density  and 
distribution  of  the  features  in  the  map  and  image.  L.N.K.  recommends  that 
some  experiments  be  done  to  estimate  the  minimum  requirements  on  the  num¬ 
ber  and  distribution  of  features. 

This  experimentation  would  proceed  in  two  ways.  One  approach  first 
involves  human  digitization  of  the  desired  features,  such  as  high  curva¬ 
ture  points  or  intersections  from  an  image.  A  map  of  the  image  would  then 
be  constructed  by  digitizing  features  from  a  different  image,  but  of  the 
same  scene,  whose  registration  with  the  first  image  is  known.  The  accu¬ 
racy  of  the  registration  process  would  then  be  examined  by  studying  the 
effect  of  varying  the  distribution  and  number  of  features.  The  experiment 
would  be  repeated  with  several  sets  of  images  so  that  more  valid  conclu¬ 
sions  may  be  reached. 

The  second  approach  is  more  elaborate.  For  an  imaginary  scene  a 
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terrain  model  would  be  constructed  wherein  every  point  would  be  repre¬ 
sented  by  its  three  coordinates  (x,y,z).  An  image  of  the  scene  would 
be  constructed  from  the  terrain  model  by  performing  a  perspective  and 
projective  transformation  on  the  map  points.  An  arbitrary  RS&T  trans¬ 
formation  on  the  image  would  then  be  calculated.  Feature  points  on  the 
map  and  image  would  be  selected  and  the  result  of  varying  their  distri¬ 
bution  and  density  on  the  registration  process  would  be  investigated. 

III.  L.N.K.  recommends  that  in  parallel  with  the  above  experiments,  a  math¬ 
ematical  analysis  of  the  registration  process  be  conducted  to  investi¬ 
gate  the  sensitivity  of  the  L.N.K.  registration  procedure  to  "noise". 

A  simplified  approximation  to  the  problem  is  to  consider  two  identical 
isolated  point  images  (such  as  high  curvature  point  images)  and  deter¬ 
mine  the  maximal  variation  in  the  ideal  rotational  and  translational 
parameters  as  a  function  of  the  size  of  a  pertubation  of  points  in  one 
of  the  images.  This  variation  will  also  depend  on  the  length  of  the  ab¬ 
stract  vectors  constructed  by  the  L.N.K.  procedure.  This  analysis  would 
be  used  in  designing  procedures  for  selecting  abstract  vectors  for  the 
L.N.K.  registration  procedure.  It  would  also  be  used  to  determine  the 
expected  "noise"  of  the  procedure,  and  the  accuracy  of  the  procedure  as 
a  function  of  its  noise. 

IV.  The  LNK  registration  procedure  requires  clustering  to  be  performed  in 
the  transformation  parameter  space.  LNK  has  only  performed  clustering 
on  three  parameters.  The  clustering  algorithm  needs  to  be  extended  to 
four  parameters  to  handle  the  full  RST  transformations.  LNK  recommends 
that  further  investigation  and  error  analysis  of  clustering  algorithms 
be  done  and  that  the  effect  of  varving  the  size  and  dispersion  of  the 
cluster  on  the  algorithms  be  investigated  experimentally  using  existing 
data. 


102 


L.N.K.  recommends  that  a  set  of  hierarchical  edge  image  registration  algo¬ 
rithms  be  devised  and  tested.  The  algorithms  would  find  successively 
smaller  windows  with  lower  resolution  from  an  original  edge  image.  Re¬ 
gistration  would  be  performed  by  matching  the  smallest  windows,  which  have 
the  lowest  resolution,  with  the  full-sized  map,  which  also  has  low  resolu¬ 
tion.  After  finding  tentative  matches  at  low  resolution  larger  windows  *.t 
high  resolution  would  be  examined  to  select  a  good  registration. 

One  method  of  matching  the  window  with  the  map  should  use  the  LNK 
registration  procedure.  A  second  method  should  use  the  sequential  correla¬ 
tion  technique.  The  hierarchical  approach  may  speed  up  the  registration 
process. 

LNK  recommends  that  interleaving  of  registration  and  feature  extraction 
be  investigated.  All  registration  procedures  presented  in  this  report 
assume  that  all  feature  extraction  has  been  done  prior  to  the  time  of 
registration.  To  reduce  the  cost  of  feature  extraction  and  registration 
it  may  be  desirable  to  determine  an  approximate  registration  from  a  lim¬ 
ited  set  of  features  and  then  use  the  approximate  registration  to  locate 
areas  suitable  for  further  feature  extraction.  These  features  can  then 
be  used  to  obtain  an  improved  registration.  This  procedure  can  be  iterated 
until  a  single  good  registration  arises. 

In  the  case  of  image  to  map  registration,  the  map  can  be  used  to 
specify  areas  in  the  image  which  are  likely  to  provide  good  edges  under 
the  approximate  registration.  For  image  to  image  registration  one  can  do 
extensive  edge  detection  in  one  image  and  use  this  to  guide  the  search  for 
good  edges  in  the  other  image.  Tne  LNK  procedure  could  be  modified  to  this 

interleaved  form. 

Several  important  questions  must  be  investigated  before  the  above 
interleaved  procedure  can  be  effectively  applied.  Methods  need  to  be 
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devised  and  investigated  for  selecting  portions  of  a  map  having  features 
which  are  likely  to  prove  useful  for  quick  registration.  This  selection 
process  will  involve  defining  various  edge  and  point  measures  such  as 
the  number  of  high  curvature  points  in  the  neighborhood  of  a  given  high 
curvature  point.  A  more  elaborate  measure  might  involve  simple  statistics 
on  the  distribution  of  abstract  vector  lengths  and  orientation  in  the 
neighborhood  of  a  point.  These  and  other  measures  could  be  used  in  the 
selection  of  the  type  of  feature  extractor  and  the  feature  extractor  para¬ 
meters  to  be  applied  in  a  given  portion  of  the  image. 
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Appendix  A  LNK  Registration  Procedure  Software 


The  LNK  registration  procedure  software,  as  existing  on  a  HP2100MX  mini¬ 
computer  at  ETL,  consists  of  one  load  module  containing  one  main  program  and 
several  subroutines.  The  program  structure  is  illustrated  below: 


COMMANDS : 


EMATCH 

i 

j  ALMTCH  j 
!  EMATCH  • 


The  registration  program  takes  a  set  of  image  edges  from  the  file  IMAGE 
and  a  set  of  map  edges  from  the  file  MAP  and  computes  the  possible  transforma¬ 
tions  that  will  transform  the  image  onto  the  map.  The  technique  used  is  to 
cluster  the  possible  transformations  in  3-dimensional  space  and  select  those 
with  the  strongest  support. 

The  data  structures  used  reside  in  the  four  common  blocks,  presented  in 
A. 1,  and  the  match  weight  matrix,  MATCHS (200,30)  is  in  the  common  block  MATRX 
which  resides  in  the  RTE-IV,  extended  memory  facility. 
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T'ne  command  to  initiate  the  i-dimensionai  clustering  that  does 
the  bulk  of  the  registration  work,  starts  in  column  one  and  is  CliUSTR. 

A  command  that  can  be  presented  independently  or  the  CL'u'STR 
command,  is  the  EMATCH  command.  This  command  allows  the  user  to  present 
a  possible  transformation  and  have  the  software  evaluate  the  strength  of 
the  transformation  given  the  edge  information  input  as  a  result  of  the 
INPUT  command  processing.  This  command  can  be  given  before,  after  or  in 
place  of  the  CLUSI'K  command.  The  format  is: 

EMATCH  <THETAl)  ,  \XS>,  <YS>  ,<DTO(> 

w‘.  1  o  r  e 

!'■  Ill i'A i  -  rotational  angle,  in  degrees,  of  the  transformation 
to  be  evaluated,  (integer) 

XS  -  X-shift.  of  transformation  to  be  evaluated  (real) 

VS  -  Y-shift  of  transformation  to  be  evaluated  (real) 

A  if)  l,  -  Jistance  tolerance  to  be  used  in  computing  strength 
of  matchs  for  this  transformation,  (integer) 

.  The  command  to  terminate  processing  is  simply  FINISH. 
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D.  1  Crocs  containing  corunon  blocks  used  in  registration 


C 

Rli  GIST  l'ROC***A*.KH.P  .  RFGISTRAT  ION  COPE  VI  'EB7V  K**  *  k 
C 

INTEGER  XKPOO)  ,X2<200)  Y  1  <200)  ,Y2<200)  ,UV<30>  , 

E  U2<30  )  ,VM  <30  )  ,  up  (3(0  .  N IMAGE  ,  N’MAP  ,  IXCF.N  T  ,  T.  YC.ENT  ,  MXCENT  , 
fc  MYCENI  ,RSlZF,CSTZr 
REAL  TMF.T<2 ,200  , R AD < 2 , 2 0 0  ) 

C 

r: 

f :  \  i  ,  X  2  ,  Y  J  ,  Y 2  --  A  P  R  A  Y  S  I-  0 R  l  NDP  0  ..  N  T  S  0 F  I M  A G E  I .  I  N L 8 

C  IM  ,  U2  ,  VI  .  V2  -  ARRAYS  FOR  ENDFMJTNTS  OF  MAP  LI  Ml  S 

C  THEY  -  ANGLES  FOR  IMAGE  AND  MAP 
RAD  -■  RAD  FOR  IMAGE  AND  MAP 

N IMAGE  -  *  OF  LINES  IN  IMAGE  <*  OF  ROWS  IN  MATRIX) 

0  Ni-fAP  f  OF  LINE'S  IN  MAP  <  #  OF  COLUMNS  IN  MAIR1X) 

G  I  <f ENT , IYCFNT  -  X  AND  Y  00  IMAGE  CENTER 

r:  MXCENT  , MYCFNT  --  X  AND  Y  OF  MAP  CENTER 

0  RSTcE.  -  *  OF  ROWS  IN  MAI  COS 

C  CG I  EE  --  #  OF  COLUMNS  IN  HATCHS 

C 

C  NAMED  COMMON  REGIS  IS  USED  TO  HOLD  J  NFC) 

C 

COMMON  /PI  GTS/  XI.  ,  X  2  .  Y  j.  ,  Y  2  ,U  V  ,Ur>  AM.  ,V2  ,  THE  I  .PAD  . 

*  N  IMAGE  ,  NMAP  ,  I XCFNT  ,  ]  YC.EN  I  AY  CENT  ,  MYCFN  I  ,  R  S  .1  EC  .C'"»  I  <!L 
r 

(;***  END  OF  M  A  C  k  0  R 1 1 G I S  T  ,  H  .PA'  L  R  S I U  N  9  F  £  B7V***** 

o 

C 

C**  **YPPCC  DUCKET  H.P.  VERSION  Vf-ElCW****** 
r 

0  DEFINES  CLUSTER  MATE  ill'  0  .  C.  US  I'ER  )  NG  IN  c-P  <  THE  TA  ,  X  .  Y  ) 
PKANT  --  CONI  AIMS  COUNT  <’i  H  OF  HITS  IN  EACH  BUCKET 
0  Pki'f.rF  -  CO  N?  A  INF  CCUPM  01  S'  OF  HITS  IN  Oil  GOT  MATRIX 

C 

INTI  I-ER  BK  TON  1(10.1  (J  OF  B!<  1  OFF  (  M>  ,  M!  ,  MM 
!  N  TE  ill  R  NcEPOKiO)  ,  NUM  i  OF  ,  NIlMX  .Nljlir  .N.’EROO'  Ml) 

Common  /r  i_str /pi;  ton  i  .  ;<;<  1 1  ?  f  ,n/lkot  ,nze  roo  .  nom  t  of  ,nuh<  .numt 

c 

r>*.*  L‘;D  01  PROC  BUCKET  Mt  k 
C 


no 


C**** * P R T F L G  MACRO  H  .  P  .  V E R S  ION  9 F E B 79  +  *** * 

c 

1 N  r  E  G  E  R  P  R  T  I N  P  ,  P  R  1'  M  '•»  M  ,  P  R  T  B  U  K  ,  P  R  1  3  B  K 
r 

PPTTNF  -  PR  [NT  INPUT  PI  AG  ( RFC  T  ANGULAR  AH!)  I’ni.AI’) 

C  PRTMM.M  -  PR  [NT  NATCH  KIT  H'.t  IT  MATRIX  FLAG 
C  PRT  I’UK  -  POINT  7-0  CLUSTERING  MATRIX  FLAG  (  BUCKL  T  B  ) 

C  PRTSBK  -  PRINT  SMOOTHED  Cl  USHERING  MATRIX 

C 

['.  NAMED  COMMON  PR  TFT.  IB  USED  TO  HOLD  FLAGS 

r 

C  DM i'i 1 I i i  / 1>  •?  T FI  /  P  R  T1  Np  ,  P  R  T  MUM  ,  P  R  T L|J |<  ,  P  R  T S  EK 
r 

C  *  *  *  END  OF  PR  l. iC  PRIILC 
0 
c 

r  *  it  .*  *  * p  r  0  r  c;  >  |  s  I)  r  0  Fi .  P  .  0  T.  R  St  ON  v  F  E  B  7  *  *  *  *  *  *  *  *  * 

P. 

INTEGER  PX(Q)  ,  DYO.T) 

C  DX.DY  -  CHANGE'S  IN  X  AND  Y  FOR  CHAIN  CODES 

C  TIP-  NAMED  COME  ON  D'  L  T  B  HUi.DS  THESE  Prtl.Uf.S  AND  1  S  1 N  |  i  l  Al  I  L I  D  T  M  Tip 

r:  to  OOP  DATA  p  R  OCR  •  iM 

c: 

COMMON  /BELTS/  D X  ,  I) Y 
C 

G*  +  *  END  or  KR'jC  CH3DF8 

c 


Ill 


ft 

t 


D. 2  Possible  commands 


The  command  needed  to  get  things  started  is  the  INPUT  command. 
The  format  is  as  follows  (starting  in  column  1): 

INPUT  <ANGTOL) ,  <STOL>  , <DTOl)>,  <NUKLE\0, <£rTBUk)>,<PRTOFF>  ,  <J?RTMWm\ 
S.PWlDTllX  <LENCI!K^>,  <y  ;IRESh)> 

Where  the  input  fields  have  the  following  significance: 

ANC.TOL  -  angular  tolerance  (in  degrees)  between  model  edge  and  a  rotated 
image  edge.  Used  in  the  routine  ALMTC.H . 

SUM,  -  segment  tolerance  for  generating  the  line  segments  in  EDGMAP. 

h POL  -  distance  tolerance  used  in  computing  strength  of  match  between 

a  model  edge  and  a  transformed  image  edge.  Calculat ion  using 
DTOL  is  i.n  routine  EMATCil . 


NL’MLEV  - 

number 

or  cluster 

ing  levels  desired. 

Maximum  currently 

poss  Lbli 

is  5 . 

Used  as  a 

controlling  parameter 

in  REGDR. 

PRTBUK  - 

flag  fo 

r  printing 

original  clustering 

matrix.  if  value 

is  1 , 

n:a  t  r  L 

is  printed 

at  each  level. 

PRTOFF  - 

f  1  a  g  f  c 

r  pr in. ting 

offset  clustering  matrix.  If  value  is 

1. 

ma  t  r  i 

is  printed 

at  each  level. 

pkt'WM  - 

flag  fu 

r  printing 

match,  weight  matrix. 

PWIDTH  - 

number 

of  columns 

available  on  output 

line.  Can  be  eitli 

er  7  2 

or  112. 

‘.FNGilK  - 

i  f  Vo  t 

Co  i  ,  on  i  v 

edges  whose  lengths 

are  approximate ly 

the 

same  w 

1  i  he  coin;1 

a  red . 

TiiRESH  - 

min  i  mur. 

:  number  of 

1  ine  reorients  th  it  r 

lust  p  ’.ss  through  a 

,  Nncket 

,’c  Lne  iiiriiest  clustering  level,  in  order  for  that  bucket  to  be 
vrotiS  :  d e  r ■ -* .i  as  a  peak.. 


112 


Kim  I L  tno 


RE  COR 


CNYTKP 


eegmap 


Arguments 


Program 


(SUMI.F.V,  SCALES. 

soaeky  ,  scalet  ,  t.'xbnd  . 

LOY EXP ,  LOTUS!)) 


(a.1,Y1,X2,  Y2, RADIUS, 
THETA ) 


( Xi ,  Y  I.  ,X2,Y2,U1 ,V1, 

.■2  ,Y2,SrO:.,TilKTA!;. 

ST  A  li.X ,  STAILY  ,  eliF.ADX , 
Slti’.M'Y ) 


(  i'AT  i'A  1  EY  ,  iiKAI'X . 
HEADY  ,  IIIKI  A  1  .  l  i : . :  >  Nl'. 
CCA!  X.l  -iX.Si.Ai.Y,  I.!)'- , 
:tYrNL"!,o(;Ti'Nn 


Futu:  t  ion 


Reads  and  analyzes  commands  and  caiis 
appropriate  subroutine.  Sets  up  initial 
values  for  variables  contained  in  com¬ 
mon  blocks  found  in  the  procs:  REGIST, 
BUCKET,  and  PRTFLG. 

Reads  in  the  edge  information  for  the 
image  and  then  the  map  from  DISK  FILES, 
i mage  edge  information  must  reside 
j.n  the  file  IMAGE::  13.  Map  edge  infor¬ 
mation  must  reside  in  MAP::18.  Computes 
centers  of  image  and  map  windows  and  places 
thorn  in  common  block.  Using  NUMLEV , 
computes  the  three  scales  at  each,  level 
and  stores  them  in  arrays  SCALEX , SCALEY 
and  SCALET.  Computes  lower  bounds  at 
first  level  and  stores  them  in  LOXBKD(l) , 
I.OYB.ND(i)  and  LOTBNDU  ) .  Image  edge 
information  stored  in  arrays  X1,X2,Y1,Y2 
of  common  block  REGIS.  Map  edge  infor¬ 
mation  into  arrays  G1,U2,V1,V2  of  same 
common  block. 

Takes  a  directed  straight  line  segment 
beginning  at  (XI, Yl)  and  ending  at 
(X2.Y2)  and  converts  it  to  polar  coor¬ 
dinates  ( RADIUS . THETA) . 

Given  an  image  edge  defined  by  (XI, Yl) 
and  (X2,Y2),  a  nap  edge  d-.  fined  by 
i.'Cl,7l)  and  (l]2,V2),  an  angle  of  rota¬ 
tion  of  image  to  map  edge  (THETAR)  and 
a  tolerance  (STOL) :  a  line  segment  in 
•.-space  which  represents  'he  constraints 
of  the  x-shift  and  the  y-shift  in  traus- 
1  orinat ions  on  tiu.se  two  edges  is 
calculated  and  is  represented  by  (FTAIi.X, 

ST  vII.Y)  and  t. SHHADX ,  SHEADY  )  . 

creates  the  clustering  matrices  used  to 
I  i  ml  the  nn>;:  I  relevant  t  r  ins  iornit  I.  i  on  s  . 

-  iuce  two  .sets  o!  .•  1  us l e r  i n g  m.rtr  im-s 
are  u.r.'d  hv  tire  approach  vie-tuia’'  and 
offset  matr  ices')  tire  rout  trie  is  called 
twice  for  each  line  segment  it  processes. 
Tie  line  segment  is  ue’iued  hv  .TA11.X, 

TAI  I.Y)  ,  (iii'.ADX  ,  HEADY  )  and  the  angle  of 


*•*--•*  —  •*< 


Arguments 


Func  t i on 


Four  ir.e 
Name 


called  to  actually  evaluate  how  good 
each  natch  is.  Taking  the  result:;  >3 f 
EMATCli,  ALMTCti  creates  the  natch  weight 
matrix  matchs  (in  common  block  MATRX) , 
the  average  match  weight,  MTC.HWT,  the 
number  of  match  ing  rows,  N'KCHRW,  the 
number  of  matching  columns,  NMCHCL,  and  in 
which)  columns  a  match  has  been  found. 

If  LENCHK=1,  only  edges  of  approximately 
the  same  length  are  compared.  PWTDTH 
is  used  by  available  output  routine. 


E MATCH  ( S  EG l , SEG2 , GRAD , TUKTA,  Given  the  transformation  (THETA , XSH  LET , 

XSi  1 1  FT ,  Y  SH 1  FT ,  I  >10  L )  YSHIFT)  ,  the  edge  in  the  image  (stored 

in  array  SEG1) ,  the  edge  in  the  nap 
(stored  in  array  SEG2) ,  the  gradient 
direction  of  the  map  edge,  GRAD ,  the 
gradient  direction  of  the  image  edge, 
THETA,  and  the  tolerance,  DTiiL; 
computes  strength  of  match  and  returns 
a  value  between  0  and  1. 


'.Corker  routines  not  shown  Ln  program  structure. 


C|7,CK  (COPE,  PEAKS,  I  > 


.  17,  XT 


(;,vn.r/.'d,cNTSi:M> 


" R TMAT 


1  XPlFF.CSiFP.NiMAGE, 

iNAP.  PMIDi'iO 


kivaif; 


s 


on 


vi 


Maintains  a  push  down  stack  that  keeps 
of  peaks  to  no  evaluated.  Each 
entrv  contains  5  fields:  THE1A  INDEX, 

X- INDEX,  Y- INDEX,  WEIGHT  and  which  matrix 
peak  appeared  in. 

Prints  either  the  original  or  offset 
clustering  matrix  depending  on  the 
value  of  C NT SUM.  Matrices  in  common 
block  BUCKET 

Prints  the  matcii  weight  matrix  created 
in  AI.MTCH  . 

Controls  t. ii a  input  and  output  to  the 
intermediate  disk  files.  Possible 
operations  arc  OPEN,  READ, WRITE  .XEVITNP  , 
CLOSE.  Can  hive  up  to  2  riles  open 
at  once. 

K, turns  1  if  I  0,  -1  if  l  0  and  0  if 
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I 


Rent  .;;o 

M.’  ™t»  Arguments  Function 


STRAIT 


)SI’K 


rich 


the  line  is  in  THETAI.  MATNifM  indicates 
whether  incrementation  is  to  be  done  in 
the  regular  matrix  (MATNL'M=  L)  or  the 
offset  matrix  (”A7NUM=2)  .  Th.e  scale  of 
the  matrix  is  defined  by  SCALX ,LOX , SCALY 
and  LOY.  The  line  is  chased  from  tail 
to  head  and  each  bucket  of  the  matrix  that 
it  passes  through  is  incremented.  Infor¬ 
mation  about  each  lin<:>  is  written  to  an 
intermediate  output  file.  Clustering 
matrices  BKTCMT  o  BKTOFF  are  in  the 
common  block  in  proc  BUCKET. 


(AX, AY , BX, BY , 

N’NELCB , STORAG.N ,  SLINKS  , 
I  FLAG) 


Given  a  line  defined  by  (AX, AY)  and 
(BX,BY)  and  whether  it  is  r.o  use  8- 
Jirectional  or  4-direction  Freeman  codes, 
STRAIT  computes  the  links  along  that 
straight  line.  The  links  are  stored  in 
the  array  STORAG(N) .  If  the  number  of 
links,  NL1NK,  is  greater  than  N,  I  FLAG 
is  set  to  indicate  an  error. 


1 FX , FY , TX , TY , NNE I Gh 
Li  ,i.:,Xl,Nd) 


Given  a  line  from  (FX,FYJ  to  (TX.TY); 
computes  the  straightest  path,  returning 
of  LI  links  in  N1  and  it  of  1.2  links 
in  N2. 


(PEAKS,  MAXPKS  ,  NL'MPKS ,  Scans  the  two  clustering  matrices, 

" :: l FI.O , PWiDTU . THRESH)  BKTCNT  and  BKTOFF,  to  find  the  MAXPKS 

highest  values.  The  indices  of  the  high 
valued  buckets  are  stored  in  the  array 
PEAKS  with  the  number  of  peaks  actually 
chosen  being  in  NT  MPKS .  THRESH  deti.n?s 
the  minimum  value  to  be  considered  as  a 
possible  peak.  D RTF LG  and  PVIPTH  provide 
information  to  tb.e  print  routines. 


>  HE  i/v ,  XSi!  1 1  l  ,  i  S d  i FT . 
ANCTOL  ,  D  VOL ,  PV.'I'JTH  , 
LFNCIil  .  MTi'i!W  P  . 

N'MCHK.i ,  NMGlh  L.COLMCH) 


Given  a  transformation  defined  by 
Ci'HETA ,  XSH  !  FT ,  YSii  1  FT)  ,  in  angle  Ldov '  nice 
in  ANGTOL  and  a  iiisi.vc.ee  tolerance  of 
P'TOL;  an  evaluation  of  how  good  the  trans¬ 
formation  is  made.  This  is  done  by 
transforming  any  image  c  dge  that  possess 
approximately  the  correct  gradient  onto 
a  map  edge.  The  Inaction  EMATCH  is 
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